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INTRODUCTION 


The  subject  of  this  research  is  image  processing  of  digitized  mammograms.  Its  purpose  is 
to  improve  the  extraction  of  information  from  mammograms  through  image  processing 
techniques  so  as  to  facilitate  better  detection  and  delineation  of  densities  and  lesions.  The  focus 
of  this  research  is  on  mammographic  density  quantification  and  lesion  delineation.  Computer- 
assisted  analysis  of  mammographic  density  would  provide  an  objective,  quantitative  measure  of 
cancer  risk  factor.  This  measure  will  be  useful  in  total  risk  analysis  in  several  ways;  in  choosing 
alternative  screening  paradigms  in  determining  the  risk-benefit-ratio  for  the  administration  of 
toxic  preventive  measures,  in  signaling  double  reading  of  the  mammograms,  and  in  selecting  the 
appropriate  lesion  detection  algorithm. 

This  research  has  the  following  main  aims. 

1.  To  develop  and  implement  a  fuzzy  object  definition  method  for  the  detection 
and  delineation  of  parenchymal  density,  masses  and  microcalcifications  in 
digitized  mammograms. 

2.  To  develop  and  implement  a  fuzzy  object  definition  method  for  the 
classification  of  lesions  and  mammographic  densities. 

3.  To  conduct  evaluation  studies  using  histologically  verified  mammographic  data 
to  determine  the  efficacy  of  the  proposed  methods  of  lesion  detection  and 
quantitative  classifications. 
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Tasks  1, 2, 3,4:  Collect  Image  Data,  Modify  Affinity  for  Vector  Valued  Features,  Implement 
and  Evaluate  Lesion  Detection  Methods. 

We  have  digitized  about  200  existing  mammograms  (at  a  resolution  of  100  microns)  from 
our  patient  database  in  the  hospital.  These  images  were  transmitted  to  our  (Medical  Image 
Processing  Group  -MIPG)  facility,  converted  to  the  format  of  3DVIEWNIX  [1]  (the  software 
platform  used  in  the  project),  and  stored  on  a  medium.  These  data  contain  normal  studies  as  well 
as  studies  with  masses  and  microcalcifications. 

We  spent  a  considerable  amount  of  time  investigating  the  affinity  relations  that  are 
appropriate  for  segmenting  digitized  mammograms.  A  novel  scale-based  fuzzy  affinity  and 
connectedness  method  has  been  developed,  implemented  in  3DVIEWNIX  and  tested.  A  brief 
summary  of  its  approach  [2,  3]  is  given  below. 

Fuzzy  affinity  between  two  nearby  pixels  is  a  reflexive  and  symmetric  relation  whose 
strength  lies  between  0  and  1  and  indicates  how  the  pixels  locally  “hang  together”  in  the  image. 
The  notion  of  fuzzy  affinity  between  two  pixels  is  expressed  as  a  non-decreasing  function  of 
fuzzy  adjacency  (dependent  only  on  the  distance  between  them),  their  homogeneity,  and  their 
agreement  to  some  global  intensity-based  object  property  or  feature.  In  determining  the 
homogeneity  and  feature-based  components  of  affinity  between  two  pixels  p\  and  P2,  a 
neighborhood  around  both  p\  and  P2  are  considered.  The  size  of  this  neighborhood,  called  scale,  is 
not  fixed  but  depends  on  the  size  of  the  largest  homogeneous  region  locally.  The  scale  is  first 
computed  automatically  for  all  pixels  in  the  image. 

Fuzzy  connectedness  between  any  two  (not  necessarily  nearby)  pixels  p\  and  pi  is  a  fuzzy 

relation  whose  strength  lies  between  0  and  1.  It  is  determined  by  considering  an  possible  paths 
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between  p\  and  pi.  A  path  is  simply  a  sequence  of  nearby  pixels.  A  strength  of  connectedness  is 
assigned  to  each  path  which  is  simply  the  smallest  affinity  of  successive  pixels  along  the  path. 
The  strength  of  connectedness  between  p\  and  pi  is  the  largest  of  the  strengths  of  all  possible 
paths  between  p\  and  pi. 

In  one  way,  we  can  imagine  that  two  pixels  are  strongly  connected  to  each  other  if  there  is  a 
path  between  them  through  locally  homogeneous  regions.  We  call  it  homogeneity-based 
connectedness.  In  another  way,  the  connectedness  between  two  pixels  could  be  imagined  as  their 
likeliness  to  fit  the  object  features  and  unlikeliness  to  fit  the  background  features.  We  call  it 
feature-based  connectedness.  One  may  note  that  feature-based  connectedness  does  not  have  the 
notion  of  path  homogeneity  and  thus  fails  to  overcome  the  effects  of  slowly  varying  components 
over  the  image.  On  the  other  hand,  homogeneity  based  connectedness  suffers  from  the  fact  that,  in 
some  applications,  objects  merge  with  background  so  smoothly  that  always  there  is  a  path  from 
object  to  background  through  locally  homogeneous  regions.  This  is  more  pronounced  when 
blurring  or  partial  voluming  is  high.  Moreover,  especially  in  a  thin  branch  of  an  object,  often, 
there  are  small  regions  with  high  inhomogeneity  that  stop  the  homogeneity-based  paths  for  the 
rest  of  the  branch.  Therefore,  the  two  different  notions  of  connectedness  are  combined. 
Additionally,  the  notion  of  scale  allows  the  correct  estimation  of  homogeneity  and  object  features 
that  is  insensitive  to  noise. 

A  detailed  mathematical  formulation  of  scale-based  fuzzy  affinity,  connectedness,  the 
associated  algorithms,  and  their  validation  on  both  2D  and  3D  clinical  images  and  phantoms  is 
presented  in  [2,  3].  A  careful  statistical  evaluation  indicates  that  the  scale-based  method  is  much 
superior  to  the  original  method  described  in  [4]. 
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We  have  investigated  several  methods  for  the  delineation  of  lesions  in  digitized 
mammograms.  We  have  basically  pursued  two  types  of  approaches. 

The  first  approach  is  based  on  fuzzy  connectedness.  It  looks  for  abnormality  in  the  network 
by  high-strength-of-connectedness  paths  within  the  image.  The  strength  of  connectedness  of 
every  connecting  path  between  every  pair  of  pixels  is  determined  using  the  method  described  in 
[3].  This  method  needs  some  further  work.  There  are  many  avenues  here  which  we  did  not  realize 
earlier.  This  approach  seems  to  offer,  without  having  to  explicitly  detect  and  delineate  lesions,  a 
method  to  identify  architectural  distortions.  One  exciting  possibility  is  to  determine  if  sufficient 
distortion  in  architecture  can  be  detected  well  ahead  of  the  time  of  appearance  of  visible  lesions. 
We  are  investigating  this  avenue  currently  and  possibly  pursue  other  funding  sources  in  the  future 
to  support  this  activity. 

The  second  approach  we  have  developed  is  called  live  wire  [5].  In  this  method,  an  operator 
initially  selects  a  point  in  the  vicinity  of  a  lesion  boundary.  At  this  time  a  “live-wire”  is  displayed 
in  real  time  as  the  operator  moves  the  mouse  cursor.  The  live  wire  represents  the  best  path  from 
the  initial  point  selected  by  the  operator  to  the  current  cursor  position.  Since  the  best  path  is 
always  computed  and  displayed  in  real  time,  the  user  can  test  how  to  select  a  largest  possible 
boundary  segment  by  moving  the  cursor  close  to  the  boundary  and  checking  how  well  the  live 
wire  snaps  onto  the  boundary.  If  this  boundary  segment  is  acceptable,  the  user  deposits  the  cursor 
which  now  becomes  the  new  initial  point  and  the  process  continues.  Typically  2-3  points  selected 
on  the  boundary  in  this  fashion  are  adequate  to  segment  and  entire  boundary.  A  3D  version  of  this 
method  has  also  been  developed  [6]  for  segmentation  of  lesions  in  MR  images.  This  method  has 
also  been  utilized  in  several  other  applications  [7,  8]. 
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Tasks  5, 6, 7:  Prepare  Reports,  Implement  and  Evaluate  Density  Classification  Methods 

A  method  for  automatically  segmenting  dense  regions  and  quantifying  density  has  been 
fully  developed,  implemented  and  tested  on  over  100  mammograms  [9,  10].  This  method  is 
described  below  in  some  detail. 

Segmentation  of  breast  from  background:  At  the  very  beginning,  using  3DVIEWNIX 
supported  LIVE-WIRE  [6]  tool,  regions  corresponding  to  pectoral  muscles  are  interactively 
excluded  when  those  are  projected  in  the  image.  In  the  entire  process,  this  is  the  only  step 
requiring  operator  intervention.  Fuzzy  connectivity  is  used  as  the  underlying  technique  in 
segmenting  breast  from  background.  To  apply  the  fuzzy  connectivity  model,  we  need  to  estimate 
several  parameters.  Studying  120  images  from  60  patients,  we  found  that  the  intensity  histogram 
always  contains  a  highly  prominent  peak  at  the  lower  intensities,  and  the  peak  is  contributed 
mostly  by  background.  The  first  prominent  peak  in  the  intensity  histogram  is  detected  and  used  to 
(roughly)  calculate  the  mean  and  standard  deviation  of  background  intensity.  To  apply  the  fuzzy 
connectivity  algorithm,  we  need  to  select  a  set  of  reference  (seed)  pixels.  For  this  purpose,  we 
assume  that  the  rightmost  column  in  the  image  always  lies  in  the  background.  We  include  all 
these  pixels  in  the  reference  set.  Fuzzy  connectedness  processing  starting  from  these  pixels  gives 
us  a  fuzzy  connectivity  image  for  background.  We  discard  connectivity  strengths  in  the  upper  half 
and  keep  the  lower  half  as  the  breast  region. 

Fuzzy  connectivity  image  for  glandular  tissue:  The  fuzzy  connectivity  method  is  used  to 
enhance  glandular  dense  regions  and  to  suppress  fat  tissues;  the  resulting  fuzzy  connectivity 
image,  in  turn,  is  used  for  automatic  segmentation  of  the  glandular  region.  The  major  task  in 
applying  the  fuzzy  connectivity  model  is  to  estimate  the  parameters  of  the  affinity  relation  and  to 
select  a  set  of  reference  pixels.  After  ignoring  the  upper  0.01  percentile  intensities,  the  mean  and 
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standard  deviation  parameters  (for  the  homogeneity  and  feature-based  affinity)  are  estimated  from 
those  parts  of  the  breast  region  falling  in  the  upper  25%  of  the  intensity  range.  Finally,  the  pixels 
in  the  breast  region  falling  in  the  upper  15%  of  the  intensity  range  are  selected  as  reference  pixels. 
The  fuzzy  connected  image  for  the  glandular  tissue  is  then  computed. 

Automatic  threshold  selection:  For  any  threshold,  the  image  is  divided  into  two  regions. 
Local  homogeneity  based  affinity  between  every  pair  of  spatially  adjacent  pixels  is  used  to  define 
their  likeliness  of  belonging  to  the  same  object  or  of  not  belonging  to  the  same  object.  The 
optimum  threshold  is  determined  from  the  associated  statistics  called  threshold  energy.  The 
threshold  with  the  minimum  threshold  energy  is  selected  as  the  optimum  threshold.  We  generate 
several  descriptors  from  the  segmented  binary  image  and  the  original  image  to  quantify  the 
glandular  tissues  as  described  below.  Note  that  all  steps  are  completely  automatic  except  the 
exclusion  of  pectoral  muscles  if  they  are  included  in  the  mammogram. 

Density  quantification:  The  method  has  been  tested  on  over  80  studies  (each  study 
produces  two  digitized  mammograms)  from  routine  exams  from  two  projections  (CC  and  MLO). 
The  population  included  normal  as  well  as  cancer  cases  (masses  and  calcifications).  Except  for 
the  exclusion  of  pectoral  muscles,  the  entire  method  has  worked  automatically  on  all  images 
wherein  all  parameters  required  by  the  algorithms  are  selected  automatically.  The  algorithms 
produced  visually  acceptable  segmentations  in  all  images.  From  the  segmented  regions  and  the 
image  intensities  in  them,  we  compute  a  set  of  density  related  parameters  including  total 
glandularity  (TG),  TG/total  fat(TF),  TG/average  fat(AvF),  TG/area  of  breast(AB),  area  of 
glandularity(AG),  AG/area  of  fat(AF),  AG/AB.  TG  and  TF  are  computed  by  integrating 
radiographic  intensity  over  respective  segmented  regions  while  AG,  AB  and  AF  are  computed  by 
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counting  the  number  of  pixels  in  the  respective  regions.  Finally,  AvF  is  computed  by  dividing  TF 
and  AF. 

To  evaluate  the  density  quantification  method,  we  tested  the  correlation  between  the 
parameters  from  the  two  projections  (CC  and  MLO).  The  correlations  for  TG,  TG/TF,  TG/AvF, 
TG/AB,  AG,  AG/AF  and  AG/AB  are  0.967,  0902,  0951,  0944,  0.959,  0.915  and  0.941 
respectively. 

We  also  conducted  a  phantom  experiment  as  follows.  A  rectangular  parallelepiped  wax 
object  was  suspended  in  a  cylindrical  water  bath  and  imaged  at  different  orientations.  The  various 
measures  based  on  integrating  intensity  in  the  segmented  object  (wax)  region  produced  more 
accurate  density  quantification  than  the  area  measures. 

Task  8,  9,  10:  Prepare  Technical  Reports/Papers,  Compare  Among  lesion  Detection 
Strategies 

The  previously  developed  density  segmentation  method  was  based  on  the  principle  of  fuzzy 
connectedness  [4],  This  method  has  been  further  extended  and  improved,  resulting  in  the 
publications  cited  in  [3,  11-13].  We  believe  that  these  theoretical  and  algorithmic  developments 
are  very  general,  constitute  a  breakthrough  in  image  segmentation,  and  have  far  wider 
applications  than  just  in  mammographic  image  processing.  Although  not  related  to  the  grant 
under  consideration,  we  have  explored  a  variety  of  applications  for  these  same  algorithms, 
including  the  clutter-free  display  of  vessels  via  MRA  [14],  Multiple  Sclerosis  lesion 
quantification  via  MRI  [15,  16],  tumor  detection  in  the  brain  via  MRI  [17],  and  the  quantification 
of  hyper-intense  lesions  in  late-life  depression  [18].  These  algorithms  have  been  compared  with 
the  earlier  versions  of  the  algorithms  and  have  been  shown  to  be  more  robust  under  noise  and 
with  improved  segmentation  effectiveness  [11-13]. 
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A  new  class  of  interactive  methods  for  lesion  detection  based  on  live  wire  has  also  been 
developed  [5].  These  are  user- steered  methods  and  are  needed  in  extremely  difficult  segmentation 
situations.  These  methods  have  been  compared  to  manual  methods  and  earlier  live  wire  methods 
and  shown  to  be  significantly  more  efficient  with  improved  precision. 

Tasks  11,  12:  Evaluate  the  Accuracy  of  Lesion/Density  Classification,  Write  Technical 
Reports/Papers. 

The  density  classification  method  has  been  tested  on  over  150  mammograms  [10]  and  has 
shown  high  consistency  (better  than  0.96  correlation)  between  the  two  projections  (MLO  and  CC) 
of  the  same  breasts.  This  method  is  now  being  utilized  in  another  project  for  density  classification 
of  women  undergoing  hormone  therapy. 

In  the  process  of  evaluating  lesion  classification  strategies  using  fuzzy  connectedness 
parameters,  we  have  realized  that  it  is  more  sensible  to  classify  the  architectural  distortions 
quantified  through  fuzzy  connectedness  prior  to  the  appearance  of  visible  mammographic  lesions. 
Based  on  40  cases  of  masses,  for  which  we  had  2-3  prior  mammograms  available  prior  to  the 
appearance  of  masses,  we  have  evidence  that  the  fuzzy  connectedness  parameters  may  be  able  to 
describe  the  architectural  distortions  that  take  place  before  visible  lesions  appear.  This  exciting 
possibility  requires  considerable  further  work  to  prove  conclusively  the  preliminary  observations. 
If  proven,  this  may  have  a  significant  impact  on  the  early  diagnosis  of  breast  cancer. 

KEY  RESEARCH  ACCOMPLISHMENTS 

•  The  development  of  a  novel  method  of  defining  the  “hanging-togetherness”  of  dense  regions 
via  scale-based  fuzzy  affinity  and  connectedness  [3].  This  framework,  we  believe, 
constitutes  a  breakthrough  in  image  segmentation.  It  has  been  further  studied  by  other 
research  groups.  We  have  extended  this  framework  considerably  [11-13]  and  demonstrated 
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that  it  is  effective  also  in  several  areas  of  medical  image  analysis  [14-18]  other  than 
mammographic  image  processing. 

•  An  interactive  method  of  lesion  segmentation  using  live  wire  [5,  6]. 

•  An  automatic,  validated  method  of  mammographic  density  quantification  and  the 
development  of  a  host  of  intensity-based  parameters  that  are  more  accurate  than  the  measure 
of  the  area  of  dense  regions.  This  method  is  currently  utilized  in  another  project  in  which 
change  in  breast  density  is  assessed  as  a  result  of  hormone  therapy. 

•  A  novel  method  of  detecting  architectural  distortions  without  explicitly  delineating  lesions. 
The  method  is  being  tested  for  its  effectiveness  in  predicting  the  onset  of  lesions.  This 
interesting  preliminary  result  indicates  that  this  method  may  have  a  significant  impact  on  the 
early  detection  of  cancer  in  the  future  before  lesions  appear. 
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CONCLUSIONS 

1.  The  new  scale-based  fuzzy  connectedness  method  is  more  robust  and  effective  than  the 
original  method.  It  is  very  effective  for  mammographic  image  segmentation  as  well  as  in 
several  other  applications. 

2.  Glandularity  is  considered  to  be  one  of  the  strongest  factors  for  breast  cancer.  Automatic 
breast  glandularity  quantification  from  digitized  mammograms  is  practical  using  the 
proposed  method.  The  method  removes  the  subjectivity  inherent  in  interactive  threshold 
selection  techniques  currently  used. 

3.  The  live  wire  method  is  effective  in  segmenting  mammographic  lesions.  It  seems  to  be  more 

robust  than  the  active  contour  methods  commonly  used.  Its  utility  is  being  evaluated  in  3D 
(MRI)  lesion  segmentation. 

4.  The  fuzzy  connectedness  method  facilitates  various  ways  of  characterizing  the  architecture 
of  the  breast.  There  are  some  preliminary  indications  that  the  distortions  in  architectures  as 
measured  by  fuzzy  connectedness  parameters  may  predict  the  occurrence  of  visible  lesions. 
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ABSTRACT 

Tangible  solutions  to  practical  image  segmentation  are  vital  to  ensure  progress  in  many  applications  of  medical 
imaging.  Toward  this  goal,  we  previously  proposed  a  theory  and  algorithms  for  fuzzy  connected  object  definition  in 
n-dimensional  images.  Their  effectiveness  has  been  demonstrated  in  several  applications  including  multiple  sclerosis 
lesion  detection/deiineation,  MR  Angiography,  and  craniofacial  imaging.  The  purpose  of  this  work  is  to  extend  the 
earlier  theory  and  algorithms  to  fuzzy  connected  object  definition  that  considers  all  relevant  objects  in  the  image 
simultaneously.  In  the  previous  theory,  delineation  of  the  final  object  from  the  fuzzy  connectivity  scene  required  the 
selection  of  a  threshold  that  specifies  the  weakest  “hanging-togetherness”  of  image  elements  relative  to  each  other 
in  the  object.  Selection  of  such  a  threshold  was  not  trivial  and  has  been  an  active  research  area.  In  the  proposed 
method  of  relative  fuzzy  connectivity,  instead  of  defining  an  object  on  its  own  based  on  the  strength  of  connectedness, 
all  co-objects  of  importance  that  are  present  in  the  image  are  also  considered  and  the  objects  are  let  to  compete 
among  themselves  in  having  image  elements  as  their  members.  In  this  competition,  every  pair  of  elements  in  the 
image  will  have  a  strength  of  connectedness  in  each  object.  The  object  in  which  this  strength  is  highest  will  claim 
membership  of  the  elements.  This  approach  to  fuzzy  object  definition  using  a  relative  strength  of  connectedness 
eliminates  the  need  for  a  threshold  of  strength  of  connectedness  that  was  part  of  the  previous  definition.  It  seems  to 
be  more  natural  since  it  relies  on  the  fact  that  an  object  gets  defined  in  an  image  by  the  presence  of  other  objects 
that  coexist  in  the  image.  All  specified  objects  are  defined  simultaneously  in  this  approach.  The  concept  of  iterative 
relative  fuzzy  connectivity  has  also  been  introduced.  Robustness  of  relative  fuzzy  objects  with  respect  to  selection 
of  reference  image  elements  has  been  established.  The  effectiveness  of  the  proposed  method  has  been  demonstrated 
using  a  patient’s  3D  contrast  enhanced  MR  angiogram  and  a  2D  phantom  scene. 

Keywords:  Image  segmentation,  fuzzy  connectivity,  object  definition,  object  delineation 

1.  INTRODUCTION 

Two-  and  higher-dimensional  images  are  currently  available  through  sensing  devices  that  operate  on  a  wide  range 
of  frequency  in  the  electromagnetic  spectrum  —  from  ultrasound  to  visible  light  to  X-  and  7-rays.1  The  activity  of 
defining  meaningful  objects  in  these  images,  generally  referred  to  as  image  segmentation,  spans  over  three  decades.2 
The  present  paper  falls  in  this  category  and  deals  with  an  extension  of  a  previous  work3  which  was  also  motivated  by 
the  problem  of  defining  objects  in  multidimensional  medical  images.  Defining  objects  in  these  image  data  is  funda¬ 
mental  to  most  image-related  applications.  It  is  obvious  that  defining  objects  is  essential  prior  to  their  visualization, 
manipulation,  and  analysis.  Even  operations  such  as  image  interpolation  and  filtering,  seemingly  unrelated  to  object 
definition,  can  be  made  more  effective  with  object  knowledge. 

Object  definition  in  images  may  be  considered  to  consist  of  mainly  two  related  tasks  -  recognition  and  delineation. 
Recognition  is  the  process  of  determining  roughly  the  whereabouts  of  the  object  in  the  image.  Delineation ,  on  the 
other  hand,  is  a  process  that  defines  the  precise  spatial  extent  and  composition  of  the  object  in  the  image.  A  variety 
of  approaches  have  been  taken  in  biomedical  imaging  applications,  wherein  the  degree  of  automation  for  recognition 
and  delineation  ranges  from  completely  manual  to  completely  automatic. 
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Among  automatic  approaches  to  recognition,  two  classes  may  be  identified:  knowledge-based  and  model-based 
Knowledge-based  methods^  form  hypotheses  relating  to  objects  and  test  them  for  recognLg  object  pTuSt 
some  preliminary  delineation  is  needed  for  forming  and  testing  hypotheses  relating  to  object  components.  Mode'i- 

based  methods  utilize  predefined  object  models  to  optimally  match  image  information  to  models  for  recognizin'- 
object  components.  recognizing 

Approves  w  delineation  may  be  broadly  classified  into  two  groups:  boundary-based  and  region-based.  Boundary- 
met  o  s  produce  a  delineation  of  the  object  boundaries  in  the  image  whereas  region-based  methods3,12-1”4 
hTS/dIme,f ti°nS  m  f0nn  0f  the  resi0n  occupied  the  °bi«*  >n  the  image.  Each  of  these  groups  mav 
tS  by  h^d  or  3  “  “d  ^  -  dePMdi,,S  “  "",e,her  the  ^  ^/boundaries  are 

alffoSh^bfeCtfmatter  °f  this  paPer  ^  elated  to  delineation.  In  a  previous  paper,3  we  described  a  theory  and 
algorithms  for  fuzzy  connected  object  definition,  treating  a  given  image  as  a  fuzzv  subset  of  the  set  of  snacial 
elements  (spels)  comprising  the  image.  The  method  is  currently  utilized  in  several  medical  imaging  applications 
including  multiple  sclerosis  lesion  segmentation  and  quantification.15'17  MR  angiography,18  and  hard  and  soft 
tissue  3D  imaging  for  craniofacial  surgery.19  6  5  y  '  411(1  sort 

In  the  present  paper,  an  extension  to  the  above  definition  of  fuzzy  objects  is  proposed.  In  the  proposed  method 
of  relative  fuzzy  connectivity,  instead  of  defining  an  object  on  its  own  based  on  the  strength  of  connectedness  all 
co-objects  of  importance  that  are  present  in  the  image  are  also  considered  and  the  objects  are  let  to  compete  amomr 
themsdvesm  having  spels  as  their  members.  In  this  competition,  every  pair  of  spels  in  the  image  will  have  a  strength 
of  connectedness  in  each  object.  The  object  in  which  this  strength  is  highest  will  claim  membership  of  the  spels  This 
approach  to  fuzzy  object  definition  using  a  relative  strength  of  connectedness  eliminates  the  need  for  a  threshold  of 
s  gth  of  connectedness  that  was  part  of  the  previous  definition.  It  seems  to  be  more  natural  since  it  relies  on  the 
fact  that  an  object  gets  defined  in  an  image  by  the  presence  of  other  objects  that  coexist  in  the  image.  AH  specified 
o  jects  are  defined  simultaneously  in  this  approach.  Its  theory  and  algorithms  are  presented  in  Sections  2  and  3 
respectively.  In  section  4,  we  introduce  the  key  ideas  behind  the  extension  of  the  relative  connectivity  model  to  an 
i  erative  amework.  In  Section  5,  we  illustrate  the  results  of  application  of  these  methods  on  a  contrast  enhanced  MR 

£23*  d?Cf  SS1'  ^  the  Same  Section-  a  2D  phantom,  we  compare  the  method  to  a  method  of  optimum 
thresholding  of  the  fuzzy  connectivity  scenes.3-0  Finally,  concluding  remarks  are  drawn  in  Section  6. 

2.  THEORY 

Terminologies  and  notations  of  our  previous  paper3  are  followed  in  this  paper.  For  completeness,  some  of  the  key 
concepts  that  are  required  in  this  paper  are  briefly  described  here.  Prior  to  this,  an  intuitive  description  of  the  key 
ideas  is  given  using  a  two  dimensional  example.  ' 

2.1.  An  Outline  of  the  Key  Ideas 

Consider  a  2D  image  composed  of  three  regions  corresponding  to  three  objects  Oi,  02,  03  as  illustrated  in  Figure 
1,  03  being  the  background.  Suppose  we  determine  an  affimv  relation3  for  each  object  that  assigns  to  every  pair 
o  near  y  spe  s  in  the  image  a  value  based  on  the  nearness  of  spels  in  space  and  in  intensity  (or  in  features  derived 
from  intensities).  Affinity  represents  local  “hanging  togetherness”  of  spels.  To  every  “path”  connecting  every  pair 
of  spels,  such  as  the  solid  curve  n  connecting  c  and  d  in  Figure  1,  a  “strength  of  connectedness”  in  every  object 
is  assigned  which  m  simply  the  smallest  pairwise  affinity  (associated  with  the  corresponding  object)  of  spels  along 
the  path  If  the  affinities  are  derived  properly,  then  in  is  likely  to  have  a  higher  strength  in  Ox  than  in  02  or  in 
3-  her,  a  path  such  as  tt2  is  likely  to  have  a  lower  strength  in  Oj  than  the  strength  of  in  Ox.  This  relative 
0  connectedness  in  different  objects  offers  a  natural  mechanism  for  partitioning  spels  into  regions  based  on 
the  strongest  paths  between  every  pair  of  spels  in  every  object.  A  spel  such  as  a  in  the  fuzzy  boundary  between  Ox 
an  3  e  claimed  by  0\  or  03  depending  on  which  pool  of  spels  it  hangs  together  more  strongly. 

2.2.  Notations  and  Definitions 

Let  T  be  any  reference  set.  A  fuzzy  subset  A  of  X  is  a  set  of  ordered  pairs  {(x,  ^(x))|x  £  X)  where  pA  :  X  [0, 11 
jthe  memiers/iip  Return  of  A  m  X.  A  fuzzy  relation  p  in  X  is  a  fuzzy  subset.  {{{x,y),pp{x,y))\x,y  £  *},  of 
XxX.  p  will  be  called  a  similitude  relation  in  X  if  it  is  reflexive  (i.e.,  Vx  £  X,pp(x,x)  =  1),  symmetric  (i.e., 


Figure  1.  Illustration  of  the  basic  ideas  of  relative  connectivity. 


Vx,y  €  X,txp(x,y)  =  and  transitive  (i.e.,  Vx,z  €  X,^(x?r)  =  maxy€^[min[/ip(x,2/),M/7(y, z)]]).  The 

analogous  concept  for  a  hard  binary  relation  is  an  equivalence  relation. 

The  pair  (Zn,a),  where  Zn  is  the  set  of  n-tuples  of  integers  and  a  is  a  fuzzy  spel  adjacency  relation  (i.e.,  any 
fuzzy  relation  that  is  reflexive  and  symmetric)  on  Zn,  will  be  referred  to  as  a  fuzzy  digital  space .  Elements  of  Zn 
will  be  called  a  spel  (short  for  spatial  element).  A  scene  over  a  fuzzy  digital  space  (Zn,  a)  is  a  pair  C  =  (C,  /)  where 
C  =  {c|  -  bj  <  Cj  <  bj  for  some  b  G  Z£},  is  the  set  of  n-tuples  of  positive  integers,  /  is  a  function  whose  domain 
is  C,  called  the  scene  domain ,  and  whose  range  is  a  set  of  numbers  [!»,#]. 

Any  fuzzy  relation  k  in  C  is  said  to  be  a  fuzzy  spel  affinity  relation  in  C  if  it  is  reflexive  and  symmetric.  In 
practice,  k  should  be  such  that  /i*(c,  d)  is  a  function  of  the  fuzzy  adjacency  between  c  and  d,  the  homogeneity  of 
their  intensities  (or  other  features)  and  their  agreement  to  some  expected  value  of  object  intensity  (or  features).20 
Further,  /x*(c,d)  may  also  depend  on  the  scale  of  the  object  at  c  and  d.20  Throughout  this  paper,  k  with  an 
appropriate  subscript  and/ or  a  superscript  will  be  used  to  denote  fuzzy  spel  affinity.  A  path  n  in  C  from  a  spel  c 
to  a  spel  d  is  a  sequence  (ci ,  C2  j  •  *  * ,  Cm)  of  m  >  2  spels  in  C,  such  that  ci  =  c  and  Cm  =  d.  The  strength  assigned 
to  a  path  \3  defined  as  the  weakest  affinity  between  successive  pairs  of  elements  along  the  path.  Thus,  the  strength 
of  7T  is  mini<i<m[^K(ci,Ci4.i)].  For  any  5  C  C,  we  say  that  the  path  t  is  contained  in  5  if  all  all  spels  in  i r  belong 
to  S.  The  strength  of  fuzzy  K-connectedness  from  c  to  d,  denoted  /i*(c,d), is  the  maximum  of  the  strengths  of  all 
paths  between  c  and  d.  We  have  shown  earlier3  that  fuzzy  ^-connectedness  is  a  similitude  relation.  Throughout  this 
paper,  the  upper  case  form  of  the  symbol  used  to  represent  a  fuzzy  spel  affinity  will  be  used  for  the  corresponding 
fuzzy  connectedness  relation. 

For  any  scene  C  =  (C,  /)  over  (Zn,  a),  for  any  fuzzy  spel  affinity  *  in  C.  and  for  any  spel  o  G  C,  the  /c- connectivity 
scene  of  o  in  C  is  the  scene  CKo  =  (C,  fKo)  such  that,  for  any  c  <=  C,  fKo{c)  =  p.K{o,  c). 


2.3.  Relative  Fuzzy  /c-Objects 
For  any  spels  o,  b  in  C.  define 

Pob«  =  {c  |  c  €  C  and  /zk(o.c)  >  mk(6,c)}.  (i) 

The  idea  here  is  that  o  and  &  are  speis  specified  in  “object”  and  “background” ,  respectively.  Note  that  Pob  =  <p  if 
6  =s  a  *  " 

A  relative  fuzzy  k -object  O  of  a  scene  C  =  (C,  /)  containing  a  spel  o  relative  to  a  background  containing  a  spei  b 
is  the  fuzzy  subset  of  C  defined  by  the  membership  function 

*f  c  ^  , 

\  0.  otherwise.  (-) 

where  tj  is  an  “objectness” function  whose  domain  is  [I.ff]  and  whose  range  is  [0, 1].  The  range  of  f(c)  is  usually  not 
[0.1]  and  /(c)  itself  may  not  directly  represent  the  degree  of  objectness.  For  example,  in  a  CT  scene  of  the  lungs, 
the  spels  consisting  of  the  interior  of  the  bronchial  tree  have  low  /(c)  values.  The  proper  choice  of  77  to  give  accurate 
values  of  the  measurements  (such  as  volume)  that  are  sought  from  the  segmented  fuzzy  object  is  a  non  trivial  issue. 
Since,  it  is  out  of  scope  of  the  present  work,  we  will  not  delve  into  this  here.  For  short,  we  will  refer  to  O  as  simply 
a  relative  fuzzy  K-object  of  C .  Some  particular  cases  are  instructive  to  study.  Suppose  /*#(&,  0)  =  1.  Then  by  (1), 
=  0  and  the  relative  /c-object  is  empty.  This  makes  sense  since  both  0  and  b  are  inside  the  “object”  in  this  case, 
and  there  is  no  meaningful  separation  between  sets  of  spels  “hanging  together”  with  0  and  with  6.  Note  also  that, 
when  C  is  a  binary  scene  and  /(o)  ^  /(&),  i.e.,  p/c(o,  b )  ^  1,  then  PQ{ is  essentially  a  connected  component  of  spels 
whose  type  is  that  of  0  that  contains  0.  To  ensure  that  this  is  a  reasonable  definition,  we  will  state  (but  not  prove) 
several  properties  of  relative  fuzzy  /c-objects.  The  most  important  among  these  is  that,  for  any  spel  p  in  Pobm  and 
most  spels  q  not  in  P0bH ,  we  get  the  same  relative  fuzzy  /c-object. 

The  following  theorem  states  the  tightness  of  relative  fuzzy  /c-objects. 

Theorem  1.  For  any  scene  C  =  (C, /)  over  (Zn, a),  for  any  fuzzy  spei  affinity  /c  in  C,  and  for  any  speis  0,  b,p  and  c 
in  C  such  that  p  €  Po6h  , 

Hk(PjC)  >  pK(b,c)  (3) 

if,  and  only  if,  c  6  PQbK- 

The  following  theorem  asserts  the  robustness  of  relative  fuzzy  /c-objects  with  respect  to  reference  spels  specified 
in  the  object  and  in  the  background. 

Theorem  2.  For  any  scene  C  =  (C,  /)  over  (Zn,  a),  for  any  fuzzy  spel  affinity  /c  in  C,  and  for  any  spels  0,  b>p  and  q 
in  C  such  that  p  €  PQb< , 

PabH  =  PpqH  If  9  €  Pbo (4) 


Note  that  the  condition  in  (4)  is  sufficient  for  P0 =  Ppq w  but  not  necessary.  The  necessary  and  sufficient 
condition  is  expressed  in  the  following  theorem. 

Theorem  3.  For  any  scene  C  =  (C,  /)  over  (Zn,  a),  for  any  fuzzy  spel  affinity  /c  in  C,  and  for  any  spels  o,  &,p  and  q 
in  C  such  that  p  €  P0bK , 

^  if,  and  only  if,  pK{b,  0)  =  /i*(g, 0).  (5) 


The  above  two  theorems  have  different  implications  in  the  practical  computation  of  relative  fuzzy  /c-objects  in 
a  given  scene  in  a  repeatable,  consistent  manner.  Although  less  specific,  and  therefore  more  restrictive,  Theorem  2 
offers  practically  a  more  relevant  guidance  than  Theorem  3  for  selecting  spels  in  the  object  and  background  so  that 
the  relative  fuzzy  /c-object  defined  is  independent  of  the  reference  spels. 

It  follows  from  Theorems  2  and  3,  by  setting  p  =  o,  that  P0qH  =  PobH  •  However,  the  constancy  of  the  relative 
fuzzy  /c-object,  even  in  this  situation  where  the  reference  spel  for  the  object  is  fixed  but  changeable  only  for  the 
background,  requires  constraints  expressed  in  (4)  and  (5). 


______  -.r  rplarivp  fii7zv  K-obieCwS.  nArneiv  that  their  domain,  that 

Thp  following  theorem  states  an  important  property  ot  relative  ruzzv  j  -  , 

me  iouowmg  tneuic  H  ^  ,  ___  twn  qnpic  0  c  c  P  ,  the  best  path  between  them  is  contained 

is  the  set  PobH ,  is  connected  in  the  sense  that  for  any  two  speis  p.  c  t  rj6, ,  h 

bv  P  b 

Theorem  4.  For  any  scene  C  =  (C.  /)  over  (Zn,  a),  for  any  fuzzy  spel  affinity  *  in  C,  and  for  any  speis  o.  b.panac 
in  C  such  that  p,  c  6  P0b. ,  the  best  path  connecting  p  and  c  is  contained  in  P0 6.  • 

3.  ALGORITHMS 

In  this  section,  we  present  an  algorithm,  named  kRFOE ,  for  extracting  the >  relative  of  a  scene 

r  =  (C  f)  containing  a  spel.  sav  o.  relative  to  a  background  containing  a  spel.  say  6.  Prior  to  this,  we  present 
another* algorithm,  nfmed \FOE.  for  creating  the  ^-connectivity  scene  of  o  in  C.  Algonthm  kFOE  is  based  on 
dynamic  programming  and  is  called  by  algorithm  kRFOE. 


Algorithm  nFOE{o) 

Input:  C  =  (C,  /),  and  k  as  defined  in  Section  2. 

Output:  x-connectivity  scene  C«o  =  ( C,/ko )•  _,1P  q  Qf  SDeis  we  refer 

Auxiliary  Data  Structures:  An  n-D  array  representmg  CKo  =  (CJko)  and  a  queue  Q  ot  speis.  We 
y  ^  fha  arrav  itself  bv  Cv*  for  the  purpose  of  the  algonthm. 


begin  .  - 

0.  set  all  elements  of  CKo  to  0  except  the  spel  o  which  is  set  to  1; 

1.  push  all  speis  c  €  C  such  that  c)  >  0  to  Q, 

2.  while  Q  is  not  empty  do 

3.  remove  a  spel  c  from  Q; 

4.  find  /max  ~  maXdgc[min[//Co(d),  d)]], 

5.  i//«K>/jf«(c)  then 

6.  set  fKo(e)  =  /m»x;  . 

7.  push  all  speis  e  such  that  pK{^e)  >  to 

endifi 

endwhile; 

8.  output  the  /c-connectivity  scene  Cx<>; 


Algorithm  0  E ( o,  b) 

Input:  C  as  (C,  /),  and  *  as  defined  in  Section  2. 

u  store  the  followms:  (1)  the  ».cohh«tmw  soene  C*.  = 
(C,  /jco),  (2)  the  ic-connectivity  scene  Cxb  =  /*&)>  and  (3)  C7- 


begin 

0.  set  CKo  =  kFOECo); 

1.  set  Cf(b  =  «irOE(6); 

2.  /or  all  c  6  C  do 

3.  if  fKo{c)  >  fKb{c)  then 

4.  set  po(c)  =*K/(c)); 

5.  else 

6.  set  /zo(c)  =  °i 
endif 

endfor 

7.  output  the  relative  fuzzy  K-object  (2; 


4.  ITERATIVE  RELATIVE  FUZZY  k-OBJECTS 
In  this  section,  we  introduce  the  key  ideas  behind  an  extension  of  relative  connectedness  to  an  iterative  framework, 
while  still  satisfying  the  key  ideas  behind  relative  connectedness  developed  in  Section  2. 


Figure  2.  Illustration  of  the  basic  ideas  of  iterative  relative  connectivity. 


Consider  the  situation  illustrated  in  Figure  2  which  demonstrates  three  objects  Ox,  0 2  and  O3.  It  is  very  likely 
that,  for  a  spei  such  as  c,  pk{o,  c)  a  /x/<-(&,  c)  because  of  the  blurring  that  takes  place  in  those  parts  where  Oi  and  O2 
come  close  together.  In  this  case,  the  stongest  path  from  6  to  c  is  likely  to  pass  through  the  “core”  of  Oi  indicated 
by  the  dotted  curve  in  the  figure.  This  core  which  is  roughly  P0b „ ,  can  be  detected  first  and  then  excluded  from 
consideration  in  a  subsequent  iteration  for  any  path  from  b  to  c  to  pass  through.  Then,  we  can  substantially  weaken 
the  strongest  path  from  b  to  c  compared  to  the  strongest  path  from  0  to  c  which  is  still  allowed  to  pass  through 
the  core.  This  leads  us  to  an  iterative  strategy  to  grow  from  o  (and  so  complementarily  from  b)  to  more  accurately 
capture  0\  (and  O2)  than  if  a  single  shot  relative  connectedness  strategy  is  used.  An  outline  of  this  formulation  is 
given  below. 

For  any  fuzzy  affinity  k  and  any  two  spels  c,  d  €  C,  define 


HKoJc,d)  =  fiK{c,d)  (6) 

P?6,  =  {c  |  c  €  C  and  nK(o, c)  >  c)}.  (7) 

Note  that  Pj*6>  is  exactly  the  same  as  P0j,,  defined  in  (1).  Assuming  that  P^1  and  x’J1  for  any  integer  i,  P‘6>  and 
Klob  are  defined  as  follows.  For  all  c,  d  €  C 


f  0,  if  c  or  d  €  P£ 

\  fiK(c,d),  otherwise, 

{c|  c€  C  and  hk{o,c)  >  tMK^(b,c)}. 


(8) 

(9) 


An  iterative  relative  fuzzy  Kl-object  O'  of  a  scene  C  =  ( C ,  /)  containing  a  spel  o  relative  to  a  background  containing 
a  spel  6  is  the  fuzzy  subset  of  C  defined  by  the  membership  function 


otherwise. 


(10) 


*  *  a„h»r,r0A  in  MR  andoeram  of  a  patient’s  left  thigh,  (a)  Whole  vessel 

via  relative  fuzzy  connectivity,  (c)  Same  as  (b)  for  vein. 

We  are  cumently  pumuing  a  rigorous  mathmnatical  and  eaperimenml  study  of  the  ponies,  robustness  end  effec- 
tiveness  of  this  Native  relative  fuzzy  connectedness  strategy  for  image  segmentatio  . 

5.  RESULTS  AND  DISCUSSION 

,  _  .  ^  r  .1  nmnr^pH  relative  fuzzv  connectedness  method  both  qualita- 

Ip  this  section,  we  demonstrate  the T®” the  resets  of  application  of  the  method  on  a  3D  contrast  enhanced 
lively  and  quantitatively.  Figm.  3  <ta*» «• the  'hole  vessel  structure  which  was  automatically  seg- 
MR  angiogram  of  a  patient  s  left  thigh.  F  gM  1  !  ,  20  „j  was  rendered  using  shell  rendering.'1 

angiography  scene  has  a  domain  of  512  x  512  x  60  and  a  voxel  size  of  0.94mm  x  0.94mm  1  80  mm. 

Figure  4  quantitatively  segmentation  of 

r^hiteSSTegion^s. carefuUy  Perfo^t^  *** 

was  then  obtained  by  assigning  to  every  pixe  in  gm  .  rtrioinal  slice  and  by  assigning  to  each  of  the  rest  of 

average  intensity  within  the  segmented  white  matter  region  m  the  original  ^ 

the  pbcels  a  constant  Intensity  equal  to  the  avmage  intensity  vWhm  th 1  segnwnted ^ ^  GauSan 

slice.  From  the  simulated  rnene,  we  crmtted  the  ■  t^Ud  comment  horn  0  to 

kernel,  (b)  a  0-mean  Gaumian  corrdated  norne  and  W  a  slmrty  mpymg  1  P  «  SCale-based  fancy 

100  across  the  columns.  Figure  4(b)  shows  the  final  tot  phantom  possible  hard 

connectivity  scene  for  the  phantom  scene  o  F  ^^h^Two  binary  scenes  is  defined  as  the 

segmentation  of  white  matter  region  was  o  .  binary  scenes  The  fuzzy  connectivity  scene 

percentage  of  normalised  £££?£?  —  id  the  white  matter 

was  then  segmented  at  a  threshold  at  which  the  ,  the  best  Dossibie  hard  segmentation  for 

scene  (i.e„  the  initial  truth)  of  Figure  4(a)  *  mmmmm  Figure  4  d  bSdSd  ^  selected  by 

3SSMT  S  reared 

ptntm  sSiTJTdoiSreoS  x  ^  %A?£dL  of  0.86  mm  *  0.86  mm.  Although,  the  final  segmentations  in 


Figure  4.  Results  of  application  of  relative  fuzzy  connectivity  on  a  2D  phantom,  (a)  Binary  white  matter  region 
manually  segmented  out  from  a  2D  slice  of  a  3D  MR  scene  of  a  multiple  sclerosis  patient’s  head,  (b)  Test  phantom 
scene  generated  from  (a)  after  adding  noise,  blurring  and  background  variation,  (c)  Scale-based  fuzzy  connectivity 
scene  for  white  matter  region,  (d)  Hard  segmented  white  matter  region  obtained  from  (c)  at  the  best  possible 
threshold,  (e)  Hard  segmented  white  matter  region  obtained  via  relative  fuzzy  connectivity, 

Figures  4(d)  and  (e)  are  visually  very  close,  the  segmentation  using  relative  connectivity  is  closer  to  the  initial  truth. 
This  observation  may  be  argued  by  the  fact  that  in  absolute  fuzzy  connectivity,  a  fixed  global  threshold  is  chosen 
while  in  relative  connectivity  the  competition  between  objects  (foreground  and  background)  is  spatially  variant. 
Effectively,  relative  connectedness  allows  a  variable  threshold  in  the  strength  of  connectivity  in  different  parts  of  the 
scene.  Most  importantly,  to  achieve  the  best  threshold  in  absolute  fuzzy  connectivity,  the  initial  phantom  truth  is 
used  which  is  not  available  in  any  real  application.  Relative  fuzzy  connectivity  eliminates  this  need. 

6,  CONCLUSION 

Based  on  our  previously  developed  framework  of  fuzzy  connectedness  and  object  definition,3,20  we  have  proposed 
an  extension  of  this  framework  that  considers  all  relevant  objects  simultaneously.  The  fundamental  premise  on 
which  this  is  developed  is  that  an  object  gets  defined  in  an  image  by  the  existence  of  other  co-objects  (including  the 
background).  We  consider  certain  regions  in  the  image  as  part  of  the  object  because  these  regions  hang-together 
more  strongly  with  object  elements  than  with  background  elements. 

One  drawback  of  the  previous  fuzzy  connectivity  theory  is  having  to  select  a  threshold  for  the  fuzzy  connectivity 
scene  to  delineate  the  object  region.  Relative  fuzzy  connectivity  provides  an  effective  and  robust  solution  to  this 
problem.  The  robustness  of  relative  fuzzy  objects  with  respect  to  reference  pixel  selection  has  been  stated.  An 
algorithm  for  computing  relative  fuzzy  connected  objects  using  dynamic  programming  has  also  been  presented.  The 
concept  of  iterative  relative  fuzzy  connectivity  has  been  introduced. 

We  have  demonstrated  the  effectiveness  of  the  proposed  method  in  artery- vein  separation  in  a  contrast  enhanced 
3D  MR  angiogram.  Based  on  a  2D  phantom  scene,  we  have  shown  that  the  segmentation  by  relative  fuzzy  connectivity 
is  better  than  that  obtainable  via  the  best  thresholding  of  the  absolute  fuzzy-connectivity  scene.  More  extensive 
experiments  in  several  ongoing  applications  axe  currently  undergoing,  as  well  as  the  development  of  th *  theory  for 
multiple  objects. 
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ABSTRACT:  The  purpose  of  this  paper  is  to  describe  a  practical 
method  for  the  clutter-free,  three-dimensional  (3D)  volume  rendering 
of  magnetic  resonance  angiographic  (MRA)  data.  In  MRA,  clutter  due 
to  artifacts  or  nearby  high-intensity  structures  prevents  clear  visual¬ 
ization  of  the  vessel  under  investigation.  We  offer  an  alternative  to  the 
manual  editing  that  is  commonly  used  to  remove  clutter.  The  method 
is  near  automatic  and  requires  the  user  to  point  at  structures  on  a  3D 
maximum  intensity  projection  (MIP)  display.  It  utilizes  recently  devel¬ 
oped  fuzzy  connected  object  delineation  algorithms  to  extract  the 
vessels  of  interest.  Because  the  resulting  definition  is  nonbinary,  it  can 
be  displayed  via  MIP  or  more  sophisticated  volume-rendering  tech¬ 
niques.  The  improved  renditions  are  illustrated  with  several  MRA 
studies.  Implementation  of  the  fuzzy  connectedness  method  proved 
to  be  effective  in  removing  the  associated  clutter  in  the  images  and, 
in  some  cases,  dramatically  improving  visibility.  Additionally,  vessels 
could  be  extracted  with  a  nominal  number  of  points  selected  within 
the  object  by  the  user  that  retained  most  of  the  information  present  in 
the  conventional  MIP  display.  This  could  all  be  performed  in  a  prac¬ 
tical  time  frame:  the  first  vessel  delineation  in  30  s,  subsequent 
delineations  in  2-10  s  per  view,  all  on  a  100-MHz  Pentium  PC. 
Automatic  delineation  of  vessels  for  3D  MRA  visualization  is  feasible 
via  fuzzy  connectedness  principles.  The  method  retains  the  original 
intensity  constitution,  an  advantage  of  the  MIP  method,  and  mostly 
eliminates  the  clutter  commonly  observed  in  MIP.  Its  speed  and 
effectiveness  make  it  feasible  for  routine  clinical  use.  ©  2000  John  Wiley 
&  Sons,  Inc.  Int  J  Imaging  Syst  Technol,  1 1 , 62-70,  2000 


I.  INTRODUCTION 

Magnetic  resonance  angiography  (MRA)  is  a  relatively  recent  de¬ 
velopment  (compared  to  X-ray  angiography).  It  has  better  sensitivity 
for  the  detection  of  vascular  abnormalities,  such  as  aneurysms,  than 
routine  magnetic  resonance  imaging  (MRI).  MRA  methods  are 
based  on  the  same  principles  as  flow  quantitation  techniques:  time- 
of-flight  (TOF)  and  phase-sensitive  principles.  Regardless  of  the 
specific  techniques  used  to  image  the  flowing  blood,  it  is  necessary 
to  isolate  the  vascular  system  from  the  surrounding  stationary  tissues 
(Stark  and  Bradley,  1992).  In  MRA,  a  major  concern  is  to  distin¬ 
guish  the  vessel  in  question  from  the  additional  clutter  introduced  by 
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artifacts.  For  the  three-dimensional  (3D)  display  of  MRA  data, 
maximum  intensity  projection  (MIP)  is  the  most  commonly  used 
technique  (Stark  and  Bradley,  1992;  Mistretta,  1993;  Owen  et  ah, 
1993;  Prince  et  ah,  1995;  Napel  et  ah,  1992;  Bluemke  and  Cham¬ 
bers,  1995).  In  MIP,  the  intensity  (shading)  assigned  to  a  pixel  in  the 
rendered  picture  represents  the  maximum  of  the  intensities  of  the 
voxels  in  the  input  volume  image  along  a  line  from  the  pixel  that  is 
orthogonal  to  the  picture  plane.  The  high-intensity  points  in  space 
are  projected  to  generate  a  MR  angiogram.  Although  this  is  concep¬ 
tually  and  computationally  a  simple  technique,  it  has  several  draw¬ 
backs.  These  include  the  presence  of  clutter  in  display,  lack  of 
information  about  structural  juxtaposition  and  constitution,  the  cre¬ 
ation  of  irregular  depth  cues,  and  the  lack  of  a  sense  of  the  real  size 
of  objects.  A  limitation  of  the  method  is  the  lack  of  selectivity 
because  every  signal  intense  structure  in  the  selected  volume  is 
depicted  (Welte  et  al.,  1996).  Nonvascular  high-intensity  patterns 
are  caused  due  to  artifacts  during  MRI  (e.g.,  susceptibility  artifacts, 
inhomogeneities  of  the  RF  fields).  Furthermore,  projections  from 
different  orientations  do  not  necessarily  contain  the  same  objects  (as 
they  may  or  may  not  be  obscured  by  another  object).  The  MIP 
procedure  along  a  projection  ray  path  is  commutative,  i.e.,  two  MIP 
images  computed  180°  apart  are  equivalent.  Each  MIP  image  thus 
has  two  equally  valid  perspectives  (correct  and  mirrored)  so  that  one 
can  only  guess  at  the  true  image  orientation  (Verdonck,  1996). 

Clutter  appears  in  MIP  displays  either  due  to  high-intensity 
artifacts,  such  as  those  stemming  from  surface  coils,  or  to  uninter¬ 
esting  vessels.  Creation  of  clutter- free  MIP  displays  has  remained  a 
challenge  in  tomographic  angiography  mainly  because  the  problem 
originates  in  the  difficult  task  of  image  segmentation.  The  successful 
solutions  to  this  problem  have  invariably  involved  some  form  of 
slice-by-slice  help  from  a  human  operator  in  defining  the  region  of 
the  vessel  in  the  given  volume  image  (Shiffman  et  al.,  1996). 

The  lack  of  information  about  juxtaposed  structures  (e.g.,  as  to 
which  vessel  is  in  front)  in  MIP  displays  is  due  to  the  lack  of  any 
structure  shading  in  the  MIP  method.  Objects  such  as  plaques  and 
calcifications  usually  have  a  graded  constitution.  MIP  cannot  capture 
this  information  in  its  displays.  Because  of  this,  various  techniques 
have  been  developed.  Such  techniques  include  TOF  and  phase- 
contrast  (PC)  MRA.  In  both  of  these  cases,  these  methods  produce 
a  projected  image  similar  to  a  conventional  angiogram.  However, 
for  analysis,  both  the  projection  and  the  individual  slices  that  make 
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up  the  image  must  be  used  (Shiffman  et  al.,  1996;  Atlas  et  al.,  1997). 
Sophisticated  rendering  strategies  that  account  for  reflection,  trans¬ 
mission,  and  emission  through  various  fuzzy  interfaces  through  the 
vasculature  are  needed  for  this  purpose. 

In  an  attempt  to  overcome  some  of  the  above  problems,  we  have 
devised  an  approach  of  vessel  definition  and  rendering  based  on  a 
recently  developed  concept  of  fuzzy  connectedness  (Udupa  and 
Samarasekera,  1996).  In  this  strategy,  an  object  (such  as  a  vascular 
tree)  is  considered  to  be  a  fuzzy  entity  where  the  object  elements 
(voxels)  “hang  together”  with  varying  degrees  of  strength.  The 
concept  and  the  associated  algorithms  have  been  effectively  applied 
quite  extensively  to  MRI  segmentation  of  the  brain  (Udupa  et  al., 
1997a,  b;  Samarasekera  et  al.,  1997;  Miki  et  al.,  1997,  1998;  Phillips 
et  al.,  1998)  and  computed  tomographic  (CT)-based  craniofacial 
applications  (Udupa  et  al.,  1997a,  b).  In  this  paper,  we  explore  their 
effectiveness  for  the  clutter-free  display  of  vessels  in  MR  A  data. 
Their  adaptation  for  MRA  is  different  from  that  for  these  other 
applications.  This  consists  of  two  steps:  fuzzy  object  definition  and 
rendering.  These  are  described  in  detail  in  the  next  section. 


II.  METHODS 

A.  Fuzzy  Connected  Vessel  Definition.  The  fuzzy  connected 
object  definition  concepts  (Udupa  and  Samarasekera,  1996)  are 
applicable  to  rc-dimensional  images.  Because  our  application  of 
these  ideas  focuses  on  3D  MRA  images,  the  description  here  will  be 
for  3D  volume  images.  We  think  of  a  3D  volume  image  /  as  a  pair 
I  =  (V,f)  where  V  is  a  3D  array  of  voxels,  and  for  any  voxel  v  in 
V,  f(v)  represents  the  image  intensity.  On  the  3D  grid  system 
defined  by  the  3D  voxel  array  V  in  a  given  3D  volume  image  /,  we 
define  a  fuzzy  adjacency  relation  a.  It  assigns  to  every  pair  (w,  v)  of 
voxels  in  V  a  degree  fxa(u ,  v)  of  adjacency  between  them.  This 
relation  a  is  intended  to  be  a  local  phenomenon  for  capturing  the 
blurring  property  of  the  imaging  device. 

Now  we  think  of  the  given  volume  image  I  as  being  defined  over 
this  fuzzy  grid  system.  We  define  over  this  fuzzy  grid  system 
another  fuzzy  relation,  called  fuzzy  affinity  k,  on  V,  that  assigns  a 
degree  of  affinity  / wK(w,  v)  to  every  pair  (w,  v)  of  voxels  in  V.  This 
relation  is  also  local,  but  it  takes  into  account  how  close  u  and  v  are 
spatially  as  well  as  in  their  image  intensities,  f(u)  and  f(v ),  and 
possibly  also  in  their  image-derived  properties  such  as  the  magni¬ 
tude  of  the  image  intensity  gradient.  For  MRA,  we  devised  the 
following  fuzzy  affinities.  These  were  arrived  at  after  extensive 
experimentation  within  our  software  framework  (Udupa  et  al., 
1994),  which  allows  combining  several  image-derived  properties 
including  the  original  intensities  to  define  affinities.  We  associate 
two  “features,”  denoted  x(u,  v)  and  c f>2(u,  v),  with  every  pair  of 
voxels  in  any  volume  image  I  =  (V,  f)  as 


4>i(u,  v)  =  ma x(f(u),f(v)). 

(1) 
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The  idea  behind  the  above  affinity  relation  is  the  following.  If  u  + 
v  and  u  and  t;do  not  share  a  voxel  face,  their  affinity  is  0.  Otherwise, 
their  affinity  depends  on  both  the  larger  of  the  two  intensities, 
denoted  (f)x(u,  v),  and  their  difference,  c />2(«,  v),  as  per  the  function 
in  Eq.  (3).  The  smaller  the  difference  and  the  greater  the  larger 
value,  the  greater  is  the  affinity.  There  are  altogether  four  parame¬ 
ters,  m , ,  crj ,  m2,  cr2 ,  that  characterize  the  affinity.  These  parameters 
are  determined  via  “training.”  The  operator  “paints”  some  typical 
vessel  regions  in  a  display  of  one  slice  of  the  volume  image  using 
the  mouse  cursor.  The  software  subsequently  determines  the  mean 
and  standard  deviation  of  the  features  <j>l  and  (f)2  within  this  region. 
The  value  of  m,  is  set  to  the  mean  for  c pl9  and  ctj  is  set  to  roughly 
four  times  the  computed  standard  deviation.  (The  logic  behind  this 
is  that  with  two  standard  deviations  on  either  side  of  the  mean,  we 
would  capture  roughly  97%  of  the  samples  in  a  Gaussian  distribu¬ 
tion.  Under  this  distribution,  therefore,  most  of  the  painted  region 
will  be  included  with  nonzero  affinity.)  The  value  for  m2  is  set  to  0, 
and  for  cr2,  four  times  the  standard  deviation  for  </>2.  Note  that  the 
training  to  learn  the  parameters  is  needed  in  only  in  one  slice.  More 
regions  may  be  painted,  but  this  is  not  necessary. 

Our  aim  is  to  define  a  3D  object,  such  as  a  vascular  tree,  as  a 
fuzzy  connected  entity.  This  is  facilitated  through  another  fuzzy 


Figure  1.  In  determining  a  fuzzy  connected  object,  the  strength  of  all 
paths  between  all  pairs  of  voxels  (such  as  u,  v)  is  considered. 
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(a) 


(b) 


(c) 


(d) 

Figure  2.  (a)  MIP  display  of  a  body  shell,  (b)  MIP  display  of  the  fuzzy  connected  vessels  extracted  from  the  body  shell  depicted  in  (a),  (c)  A 
volume  rendition  of  the  fuzzy  connected  vessels,  (d)  Interactive  slice  selection  guided  by  an  MIP  display.  The  MRI  corresponding  to  the  depicted 
slice  is  displayed. 


relation  defined  on  V ,  called  fuzzy  connectedness,  K.  This  is  a 
global  phenomenon  as  compared  to  a  and  k.  The  strength  of 
connectedness,  v),  between  any  two  voxels  u  and  v  is 

determined  as  follows.  There  are  many  possible  paths  between  u  and 
t>,  where  each  path  is  a  sequence  of  nearby  voxels  starting  from  u 
and  ending  on  v.  Figure  1  shows  three  paths  among  the  possible 
paths  between  u  and  v.  Each  path  has  a  strength  associated  with  it, 
which  is  simply  the  smallest  of  the  affinities  between  successive 
voxels  in  the  path.  Finally,  the  strength  of  connectedness  /xA (//,  v) 
between  u  and  v  is  the  largest  of  the  strengths  of  all  paths  between 
u  and  v. 

A  fuzzy  object  of  strength  x  is  a  pool  of  voxels  (together  with 
their  strength  values)  such  that,  between  any  two  voxels  u  and  v  in 
the  pool,  the  strength  of  connectedness  fxK(u,  x>)  ^  x,  whereas 
between  any  two  voxels  v  and  >r  such  that  u  is  in  the  pool  and  w  is 


not.  the  strength  fxK{u,  t/)  <  x.  Typically,  a  small  value  of  x,  no 
greater  than  0.1,  is  adequate  to  separate  the  object  from  the  rest  of 
the  image. 

For  completeness,  we  present  one  algorithm  for  fuzzy  connected 
object  tracking.  More  details  on  the  algorithm  and  other  algorithms 
and  the  analysis  of  their  behavior  can  be  found  in  Udupa  and 
Samarasekera  (19%).  This  algorithm  takes  as  input  an  MRA  volume 
image  /  =  (V,  /),  a  voxel  o  specified  in  the  vessel  region,  affinity 
k  (as  defined  in  Eqs.  1-6),  and  a  threshold  a\  Its  output  is  a  fuzzy 
object  O  represented  as  a  volume  image  O  —  (V,  /„),  such  that  for 
any  voxel  x>  in  the  pool  corresponding  to  the  fuzzy  object,  fa(v)  — 
f{v),  and  =  0  if  x>  is  not  in  the  nook  The  algorithm  uses  a 
queue  Q  and  the  volume  image  O  as  auxiliary  data  structures.  It  is 
a  slightly  modified  version  of  algorithm  kxFOE  of  (Udupa  and 
Samarasekera,  1 996). 
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ALGORITHM  kxBFT 

begin 

set  f0(c)  =  0  for  all  voxels  c  of  V; 
initialize  queue  Q  to  contain  just  o\ 
mark  o\ 
repeat 

dequeue  a  voxel  q  from  Q ; 
for  each  unmarked  6-neighbor  w  of  q  do 
if  fJbK(q ,  w)  ^  x  then 
mark  w; 

enqueue  won  <2; 
endif; 
endfor; 

until  Q  is  empty; 
for  each  marked  voxel  c  do 
set  fjc)  =  /(c); 
endfor; 
end 

B.  Fuzzy  Object  Rendering.  The  underlying  rendering  princi¬ 
ple  is  based  on  a  method  called  shell  rendering  (Udupa  and  Odhner, 
1993).  In  shell  rendering,  potentially  every  voxel  v  in  I  contributes 
to  the  rendered  picture.  This  contribution  is  influenced  in  three  ways: 
(1)  reflection  of  light  from  v  depending  on  the  strength  of  the  surface 
contained  in  v ,  (2)  emission  depending  on  the  degree  of  object 
membership  of  v,  and  (3)  attenuation  depending  on  the  opacity 
assigned  to  v.  The  shell  rendering  concept  allows  a  continuum  from 
only  reflection-based  hard  surface  rendering  to  a  variety  of  degrees 
of  translucent  renditions  that  mix  reflection,  emission,  and  transmis¬ 
sion  in  different  proportions. 

A  shell  is  a  data  structure  that  allows  rapid  rendering  of  fuzzy 
objects  by  retaining  only  those  voxels  in  its  representation  that 
make  nonnegligle  contribution  to  the  rendered  picture.  With  each 
voxel,  several  other  items  of  information  are  stored,  including  the 
voxel’s  degree  of  membership  in  the  object,  the  magnitude  of  the 
gradient  at  the  voxel  computed  from  the  given  volume  image,  and 
the  direction  of  the  gradient.  In  MRA,  for  example,  the  voxels 
outside  the  body  and  even  those  far  away  from  the  vessel  bound¬ 
ary  need  not  be  stored  in  the  shell.  However,  before  the  vessels 
are  actually  detected,  we  cannot  be  sure  as  to  which  voxels  inside 
the  body  do  not  belong  to  the  vessels.  Therefore,  we  first  create 
a  shell  that  stores  only  the  voxels  that  are  inside  the  body, 
together  with  their  descriptions.  These  voxels  are  easily  identified 
by  thresholding.  For  further  reference,  we  will  call  this  a  body 
shell.  This  shell  is  created  automatically  without  loss  of  any 
relevant  information  when  the  image  data  are  transferred  from 
the  MR  scanner  to  the  viewing  workstation  via  our  picture 
archiving  and  communication  system. 

We  then  create  an  MIP  display  of  the  body  shell,  with  3D 
orientation  selections  under  user  control.  We  utilize  this  display, 
although  cluttered,  to  guide  the  user  in  selecting  the  vessel  structures 
that  need  to  be  extracted  from  the  body  shell.  This  selection  is  done 
by  pointing  the  cursor  at  a  vessel  structure  and  clicking  the  mouse 
button.  By  this  action,  the  user  specifies  a  voxel  in  the  array  V.  This 
specification  is  made  possible  by  storing  the  coordinates  of  the 
maximum  intensity  voxel  that  contributed  to  the  rendition.  Several 
such  voxels  may  be  specified.  Typically,  a  voxel  should  be  specified 
for  each  separate  vessel  structure  that  is  either  not  connected  to  or 
loosely  connected  to  other  structures  for  which  voxels  have  already 
been  specified.  Once  this  specification  is  completed  and  the  process 


to  detect  fuzzy  connected  objects  is  initiated,  the  fuzzy  connected 
objects  containing  the  specified  voxels  are  first  detected  using  algo¬ 
rithm  kxBFT.  These  are  then  converted  to  a  shell  representation  and 
then  displayed  using  MIP. 

C.  Image  Data.  We  utilize  five  patient  MRA  data  sets  to  dem¬ 
onstrate  the  effectiveness  of  the  fuzzy  connected  object  detection 
and  rendering  methods.  These  data  sets  were  acquired  as  a  stack  of 
rapid  2D  gradient  echo  images  using  TOF  effects  to  produce  vessel 
images  brighter  than  adjacent  stationary  structures.  The  matrix  size 
was  256  X  256,  with  the  number  of  slices  varying  between  86  and 
140.  The  voxel  sizes  in  these  data  sets  ranged  from  1.02  X  1.02  X 
3  mm  to  1.17  X  1.17  X  3  mm. 

All  algorithms  were  implemented  within  the  3DVIEWNIX  soft¬ 
ware  system  (Udupa  et  al.,  1994)  and  all  experiments  were  per¬ 
formed  using  this  system. 

III.  RESULTS 

The  implementation  of  the  methods  has  been  optimized  for  routine 
interactive  use  in  the  clinical  vascular  imaging  section.  In  a  routine 
use,  the  following  steps  are  involved.  The  timings  reported  are  all 
for  a  100-MHz  Pentium  PC  with  256  MB  RAM. 

Step  1:  The  image  data  are  transferred  to  the  workstation  and 
the  body  shell  is  created  simultaneously  and  automati¬ 
cally.  This  typically  takes  1  min. 

Step  2:  MIP  display  of  the  body  shell  is  rendered.  This  takes 
5-10  s  per  view.  Figure  2(a)  gives  an  example. 

Step  3:  One  or  more  points  are  specified  on  the  vessels  in  this 
MIP  display.  A  fuzzy  connected  object  of  a  specified 
strength  x  (chosen  on  a  slider)  containing  the  points  is 
extracted  and  rendered  using  MIP  display.  This  opera¬ 
tion  takes  20  s.  Figure  2(b)  shows  the  fuzzy  connected 
vessels  extracted  from  the  body  shell  in  Figure  2(a). 

Step  4:  If  more  vessels  are  to  be  selected,  additional  points  are 
specified.  Subsequently,  the  fuzzy  connected  object 
containing  all  points  specified  so  far  is  extracted  and 
rendered  via  MIP.  This  step  takes  about  2-5  s. 

Step  5:  The  fuzzy  connected  vessel  structure  is  visualized  at  a 
fraction  of  a  second  per  view  via  MIP  or  more  sophis¬ 
ticated  volume  rendering  that  takes  into  account  trans¬ 
mission,  reflection,  and  emission.  Figure  2(c)  shows 
such  a  rendition  of  the  vessels  in  Figure  2(b). 

Step  6:  For  a  closer  scrutiny  of  selected  regions  in  the  vessel,  cut 
planes  of  arbitrary  orientation  are  interactively  selected 
and  stepped  along  vessels  to  determine  the  image  inten¬ 
sity  values  on  those  planes.  Figure  2(d)  illustrates  this 
operation.  The  operation  takes  5-10  s  the  first  time  and 
subsequently  about  1-2  s. 

Figure  3  shows  different  renditions  of  the  five  data  sets.  The  rows 
correspond  to  the  different  data  sets.  The  first  column  shows  MIP 
renditions  of  the  five  body  shells.  The  second  column  shows  MIP 
renditions  of  the  fuzzy  connected  vessels  extracted  from  the  body 
shells.  All  data  sets  required  one  to  three  points  to  extract  the  vessels 
shown. 

IV.  DISCUSSION 

As  seen  from  all  patient  study  examples,  fuzzy  connected  vessel 
definition  is  remarkably  effective  in  removing  the  clutter,  with  a 
dramatic  improvement  in  visibility  in  some  cases.  Most  of  the 
vessels  seen  in  the  MIP  rendition  of  the  body  shell  are  included  in 
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Figure  3.  Column  1 :  MIP  renditions  of  the  body  shells  of  the  five  data  sets.  Column  2:  MIP  renditions  of  the  fuzzy  connected  vessels  extracted 
from  the  body  shells.  Column  3:  Volume  rendering  of  the  fuzzy  connected  vessels. 


the  fuzzy  connected  object.  We  emphasize  that  in  the  fuzzy  con-  connectedness  >:  x)  for  the  same  study.  Here  |  •  |  denotes  the 

nected  objects,  all  original  MR  intensity  values  that  matter  (namely,  cardinality  of  the  set.  |X,  -  Xf\  represents  the  number  of  voxels  in 

those  in  the  vessels)  are  retained  in  the  fuzzy  connected  object.  In  X , that  are  not  in  Xf.  Clearly,  the  extreme  values  of  this  fraction  are 

order  to  test  the  validity  of  this  claim,  an  operator  outlined  the  major  0  and  100,  assuming  that  |X,|  ^  0  for  the  study.  The  values  of  this 

vessel  regions  in  every  slice  for  three  among  the  five  studies  (shown  fraction  for  the  three  studies  were  0.1487,  0.2201 ,  and  0.3174. 

in  columns  1,  2,  and  3  in  Fig.  3).  Denoting  the  set  of  voxels  defined  The  most  time-consuming  step  in  our  approach  after  creating  the 

in  this  fashion  for  a  study  by  X,,  we  computed  the  fraction  (|X,  -  body  shell  is  the  first  fuzzy  connected  vessel  definition  step,  requir- 

Xyj / 1 X j | )  X  100,  where  X(  denotes  the  set  of  voxels  determined  by  ing  about  20  s.  The  computations  involved  in  algorithm  kxBFT  are 

the  algorithm  to  be  in  the  fuzzy  object  (i.e.,  with  strength  of  both  computation  and  memory  intensive.  With  64  MB  RAM  (as 
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Figure  3.  ( Continued ) 


opposed  to  256  MB),  the  time  requirement  for  this  jumps  to  2  min, 
becoming  noninteractive  and  unsuitable  in  a  clinical  vascular  imag¬ 
ing  setting.  Interestingly,  on  a  300-MHz  Pentium  PC  (running  under 
Linux)  with  256  MB  RAM,  all  operations  described  in  Steps  1-6  ran 
at  interactive  speeds.  This  is  particularly  relevant  to  the  cost-con¬ 
scious  health  care  environment  considering  how  inexpensive  the 
Pentium  PCs  are  today. 

Volume-rendering  operations  that  take  into  account  transmission, 
emission,  and  reflection  allow  more  flexibility  for  portraying  the 
heterogeneity  of  MR  signal  intensities  in  the  vessels  than  the  MIP 
method.  MIP  rendering  does  not  portray  the  geometric  relationships 
among  and  the  shapes  of  the  vessels  properly  because  of  the  lack  of 
shading  coming  from  the  above  components.  This  is  demonstrated  in 
Figures  4(a)-(d),  which  show  an  MIP  rendition  of  a  fuzzy  connected 
vessel  and  its  three  volume  renditions  at  different  combinations  of 
the  three  properties.  Note  how  the  heterogeneity  of  MR  image 
intensities  comes  through  in  the  displays,  especially  in  Figure  4(d). 
The  extra  computational  time  (1-2  s)  for  such  renditions  over  MIP 
for  vessels  is  negligible. 

One  possible  disadvantage  of  fuzzy  connected  MIP  and  volume 
rendering  over  direct  MIP  rendering  of  the  body  shell  is  that  the 
contextual  information  coming  from  the  faintly  portrayed  body 
contour  in  the  latter  display  is  lacking  in  the  former.  Such  informa¬ 
tion  is  potentially  useful  in  providing  constantly  an  orientation  to  the 
viewer  for  extremely  unfamiliar  viewing  directions.  The  body  sur¬ 
face  information  is  easily  grafted  into  the  fuzzy  connected  renditions 
by  segmenting  the  body  at  the  threshold  at  which  the  body  shell  was 


created  and  by  merging  the  shells  corresponding  to  the  body  surface 
obtained  in  this  fashion  and  the  fuzzy  connected  vessels.  Figure  5(b) 
illustrates  a  volume  rendition  of  a  composite  shell  obtained  in  this 
fashion.  Figure  5(a)  shows,  for  comparison,  an  MIP  rendition  of  the 
original  body  shell. 

There  is  no  guarantee  in  our  system  that  with  a  single  seed,  all 
aspects  of  the  vessels  depicted  in  the  MIP  rendition  are  delineated  by 
the  fuzzy  connectedness  method.  To  include  vessels  that  are  left  out, 
the  user  simply  has  to  select  more  seed  points  in  such  vessels  in 
subsequent  stages  after  verifying  the  resulting  shell  renditions,  as 
described  in  Step  4. 

We  have  not  tested  the  system  specifically  for  its  ability  to  detect 
stenosis  better  than  in  the  original  MIP  renditions.  This  requires 
further  work  involving  observer  studies  and  comparison  with  X-ray 
angiography.  However,  as  indicated  by  our  experiment  comparing 
manual  tracing,  our  method  loses  less  than  0.5%  of  the  voxels 
identified  by  an  operator.  This,  combined  with  the  fact  that  our 
method  retains  all  original  intensity  information,  indicates  that  the 
MIP  renditions  of  the  delineated  vessels  contain  almost  all  the 
information  contained  in  the  original  MIP  renditions  minus  the 
clutter.  The  proof  of  this  claim  requires  further  validation. 

The  fuzzy  objects  delineated  by  our  method  may  be  utilized  for 
measuring  vessel  diameters.  Such  tools  already  exist  in  3D  VIEW- 
NIX.  However,  these  measurements  need  to  be  carefully  calibrated 
to  ensure  that  they  agree  with  truth.  We  have  not  done  this  testing  in 
this  work. 
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(a) 


(b) 


(c)  (d) 

Figure  4.  MIP  rendition  of  a  fuzzy  connected  vessel  (a)  and  its  volume  renditions  (b-d)  at  different  high  settings  of  transmission  (b),  emission 
(c),  and  reflection  (d). 


V.  CONCLUDING  REMARKS  utilizes  information  from  all  connecting  paths  between  all  possible 

We  have  presented  an  approach  and  a  system  for  removing  clutter  in  pairs  of  voxels.  We  have  demonstrated  that  the  approach  is  practi- 

MRA  with  minimal  user  effort  as  an  alternative  to  the  slice-by-slicc  cally  viable,  requiring  less  than  15  s  for  all  operations  and  less  than 

removal  of  obscuration  that  is  currently  practiced.  The  approach  is  4  s  for  all  rendering  operations  on  a  300-MHz  Pentium  PC.  We  have 

based  on  a  theory  of  fuzzy  connectedness  of  object  regions  that  not  carried  out  rigorous  evaluation  studies  for  specific  clinical  tasks 
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(a)  (b) 

Figure  5.  (a)  MIP  rendition  of  the  body  shell  of  an  MRA  data  set.  (b)  A  combined  display  via  volume  rendering  of  the  body  surface  and  fuzzy 
connected  vessels  extracted  from  the  body  shell  in  (a). 


comparing  the  MIP  method  commonly  used  in  clinical  MRA  today 
with  the  proposed  fuzzy  connected  strategies.  However,  we  have 
presented  evidence  based  on  five  clinical  studies  that  the  proposed 
method,  at  the  least,  retains  most  of  the  information  present  in  the 
conventional  MIP  display,  removes  most  of  the  obscuring  clutter, 
and  possibly  portrays  this  information  better  than  conventional  MIP. 
Although  we  have  not  experimented  with  CT  angiography,  we 
believe  that  the  approach  presented  here  is  applicable  to  CT  angiog¬ 
raphy  for  separating  vessels  from  high-intensity  structures,  such  as 
bones.  We  are  currently  investigating  ways  to  separate  arteries  and 
veins  in  MRA  images  using  fuzzy  connectedness  principles. 
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abstract 

Multiple  Sclerosis  (MS)  is  an  acquired  disease  of  the  central  nervous  system.  Subjective  cognitive  and  ambulatory  test  scores 
on  a  s'**1*  called  EDSS  are  currently  utilized  to  assess  the  disease  severity.  Various  MRI  protocols  are  being  investigated  to 
study  the  disease  based  on  how  it  manifests  itself  in  the  images.  In  an  attempt  to  eventually  replace  EDSS  by  an  objective 
measure  to  assess  the  natural  course  of  the  disease  and  its  response  to  therapy,  we  have  developed  image  segmentanon 
methods  Hawf  on  fuzzy  connectedness  to  quantify  various  objects  in  multiprotocol  MRI.  These  include  the  macroscopic 
objects  such  as  lesions,  the  grey  matter  (GM),  white  matter  (WM),  cerebrospinal  fluid  (CSF),  and  tain  parenchyma  as  well 
as  the  microscopic  aspects  of  the  diseased  WM.  Over  1000  studies  have  been  processed  to  date.  By  far  the  strongest 
correlations  with  the  clinical  measures  were  demonstrated  by  the  Magnetization  Transfer  Ratio  (MIR)  histogram  parameters 
obtained  for  the  various  segmented  tissue  regions  emphasizing  the  importance  of  considering  the  microscopic/diffused  nature 
of  the  disease  in  the  individual  tissue  regions.  Brain  parenchymal  volume  also  demonstrated  a  strong  correlation  with  the 
cBnifai  measures  indicating  that  brain  atrophy  is  an  important  indicator  of  the  disease.  Fuzzy  connectedness  is  a  viable 
segmentation  method  for  studying  MS. 

Keywords:  Image  segmentation.  Multiple  Sclerosis,  MR  imaging,  fuzzy  connectedness. 


1.  INTRODUCTION 

MS  is  an  acquired  disease  of  the  central  nervous  system  (CNS)  characterized  primarily  by  multifocal  inflammation  and 
destruction  of  myelin.  Inflammation  and  edema  are  accompanied  by  different  degrees  of  demyelin aticn  and  destruction  of 
oligodendrocytes,  and  may  be  followed  by  remyelinarion,  axonal  loss'  and/or  gliosis.  The  disease  was  first  characterized  by 
Charcot1  and  has  since  been  investigated  extensively.  The  highest  frequency  of  MS  occurs  in  northern  and  central  Europe, 
Canaria,  the  USA,  and  New  Zealand  and  South  of  Australia.2  In  the  US,  it  affects  approximately  350IXX)  adults  and  stands 
out  as  the  most  frequent  cause  of  non-traumatic  neurologic  disability  in  young  and  middle-aged  adults.  In  its  advanced  state, 
the  riiwasp;  may  severely  impair  the  ambulatory  ability  and  may  even  cause  death. 

MS  is  usually  classified  into  three  subtypes:  (1)  Relapsing-remitting  (RR):  Clearly  defined  disease  relapses  with  full 
recovery,  or  with  sequelae  and  residual  deficit  upon  recovery,  but  there  is  no  disease  progression  between  relapses.  (2) 
Secondary-progressive  (SP):  Initial  RR  disease  course  followed  by  gradual  progression  with  or  without  occasional  super¬ 
imposed  relapses.  (3)  Chronic-progressive  (CP):  Gradual  progression  from  the  onset  with  occasional  plateaus  and  temporary 
minor  improvements.  The  most  commonly  used  scale  to  clinically  measure  disease  progression  in  MS  is  the  Expanded 
Disability  Status  Scale  (EDSS)  introduced  by  Kurtzke.4  This  scale  extends  from  0  to  10,  and  is  based  on  the  functional 
systems  (visual,  brainstem,  pyramidal,  cerebella,  sensory,  bowel/bladder,  and  cerebral)  in  the  lower  scores,  and  on 
ambulation  in  the  higher  scores.  The  clinical  quantification  of  disease  severity  is  subjective  and  sometimes  equivocal.  The 
development  of  new  treatments  demands  objective  outcome  measures  for  relatively  short  trials.  Therefore,  MR  imaging  has 
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become  one  of  the  most  important  paraclinical  tests  for  diagnosing  MS  and  in  monitoring  disease  progression  in  MS.  As 
described  in  the  rest  of  this  section,  a  variety  of  MRI  protocols  have  been  utilized  in  understanding  this  disease. 

1.1.  T2- weigh  ted  imaging  (T2WI) 

The  sensitivity  of  T2WI  in  the  detection  of  MS  has  long  been  recognized.1-  6  MS  lesions  appear  hyperintense  in  T2WI.  The 
clinical  course  of  MS  is  mainly  defined  by  the  "baseline"  disability  of  the  patient,  also  called  the  disease  "burden".  The  extent 
of  lesions  on  T2WI,  expressed  usually  by  the  total  volume  of  lesions  (T2LV),  is  currently  regarded  as  the  MRI  measure  of 
this  disease  burden.  Usually  dual  echo  T2  and  proton  density  (PD)  images  are  acquired  and  the  lesions  are  detected  using  the 
combined  information  presented  in  the  two  images.  The  positive  effects  of  three  currently  available  injectable  medications 
(Betaseron,  Avonex,  Copaxone)  on  T2  lesion  load  (T2LV)  have  been  demonstrated.7 "°  Despite  impressive  progress  in  the 
treatment  and  understanding  of  MS,  both  its  etiology  and  pathogenesis  still  remain  incompletely  understood.  An  important 
limitation  of  T2WI  is  its  low  pathological  specificity.11  All  histopathological  features  of  MS,  such  as  inflammation,  edema, 
demyelination,  axonal  loss,  gliosis  and  remyelination  are  seen  as  hyperintense  lesions  on  T2WI  while  not  all  of  these 
processes  are  clinically  relevant  Reports  from  several  research  groups  have  shown  only  weak  correlations  between  clinical 
measures  of  disability  and  the  traditional  T2  lesion  load.12'14  There  are  also  contradictory  results  among  these  reports  -  for 
example,  that  T2LV  increases  consistently  and  significantly  in  RR  MS13- 13  and  that  it  does  not  ^ 16  The  reasons  for  such 
disagreements  include  imaging  and  processing  methodological  differences,  subjectivity  of  EDSS  and  the  non-accountability 
of  the  microscopic  disease,  hi  short,  there  is  increasing  evidence  that  T2LV  is  not  a  reliable  measure  of  disease/disability  in 
spite  of  its  use  in  clinical  trials. 

1.2.  Tl-weighted  imaging  with  gadolinium  enhancement  (TTWIE) 

In  MS,  gadolinium  (-DTPA)  enhancement  on  MR  images  is  characterized  by  blood  bain  barrier  (BBB)  breakdown  and 
intense  inflammation  which  represents  the  acute  stage  of  the  evolution  of  MS.17, 11  hi  T1WEE,  gadolinium  enhancement 
appears  as  a  homogenous  strongly  hyperin  tense  lesion  or  as  a  ring-shaped  hyper-intensity  at  the  edge  of  chronic  reactivated 
lesions.  Enhancement  is  more  sensitive  than  either  clinical  examination  or  T2WI  in  detecting  disease  activity  and  potentially 
can  separate  clinical  groups.19  Delayed  imaging  (50-60  minutes  after  gadolinium  injection),  triple  doses  of  pdolinium,  and  a 
single  dose  with  MT  saturation  to  suppress  normal  brain  have  all  been  shown  to  increase  the  number  of  detectable  enhancing 
lesions.  (There  are  disadvantages  to  high-dose  contrast  studies  including  possibly  more  false  lesions,  more  scan  time  and 
more  expense.)  These  results  suggest  that  enhancement  is  not  an  "all  or  none"  phenomenon.  Enhancement  may  precede  the 
development  of  T2  hyper-intensity  and  clinical  symptoms,  which  suggest  that  the  BBB  abnormality  may  be  the  crucial  event 
in  the  inflammatory  cascade.13,21  Further,  we  know  that  there  is  microscopic  manifestation  of  the  disease  in  normal  appearing 
white  matter  (NAWM)  which  is  not  apparent  under  visual  examination  of  the  images. 

1J.  Magnetization  transfer  imaging  (MTI) 

MS  is  now  believed  to  be  a  diffuse  process  that  extensively  involves  the  white  matter  and  is  not  restricted  to  the  focal  regions 
of  disease  activity  visible  as  "lesions"  in  conventional  T2,  PD  and  T1WIE  images.  This  notion  came  about  mainly  from 
magnetization  transfer  (MT)  image  analysis,22  but  subsequently  also  supported  by  postmortem23  revealing  microscopic 
disease  characterized  by  edema,  cellular  infiltration  and  demyelination  in  macroscopically  normal  appearing  white  matter.  In 
MT  imaging,24  two  consecutive  sets  of  images  are  required,  one  with  off-resonance  pre-saturation  of  the  relatively  immobile 
macromolecular  protons  and  one  without.  A  magnetization  transfer  ratio  (MTR)  image  is  then  computed  from  the  two  images 
according  to  (Ma  -  A/J  /  where  M,  and  M0  are  pixel  intensities  with  and  without  pre-saturation.  Published  studies22'27 
suggest  that  the  MTR  values  for  white  matter  in  normal  subjects  are  highly  reproducible  within  each  institution  with  a 
variation  of  less  than  2%.  Several  studies  have  demonstrated  that  the  microscopic  aspect  of  the  disease  can  be  characterized 
by  the  MTR  values  measured  in  different  regions  (WM,  lesions,  the  periphery  of  lesions)  in  a  much  more  specific  fashion 
than  by  T2LV. 

1.4.  Tl-weighted  imaging  (T1WI) 

Another  standard  MRI  technique  that  is  routinely  used  in  MS  is  T1WI  (SE).  It  often  shows  hypo-intense  areas  which  were 
suggested  to  represent  areas  of  axonal  loss  and  gliosis.23  In  a  recent  quantitative  study,29  a  significant  correlation  was 
reported  between  increase  in  disability  and  increase  in  the  volume  of  these  lesions  (T1LV).  In  another  study  involving  a 
larger  cohort  of  patients30  including  RR  and  SP  patients,  a  strong  correlation  was  found  oniy  in  the  SP  group  between 
increase  in  T1LV  and  increase  in  disability  over  3  years  of  follow-up  that  was  higher  than  the  increase  in  T1LV.  In  RR 


patients. many  T2  change  are  tocompanied^ 

^ISSSra  letoToad  to  clinica.  disabiiity.  The  obsenmd  differences  between  RR  and  SR 
groups  probably  relate  to  a  quantitative  difference  in  repair  mechanisms. 

have  been  devefooinz  MR  image  analysis  methods  for  MS  since  1993.  Our  approaches  are  guided  by  two  premises.  First, 

•  reeions  in  any  anatomic  object  have  a  heterogeneous  composition  with  or  without  the  disease.  Treating  them  m  a 
JJJL  fiLi™  bv  considering  a  voxel  to  contain  either  100%  or  0%  tissue  material  is  unrealistic  and  inaccurate.  The 
^STmediSs  fflTthe  quantitative  measures  derived  from  them  that  we  have  devised  take  into  account  the 
ES  »?^terin«s  inherent  in  tissue  regions.  Second,  as  we  have  seen  from  the  description  in  this  section 
different  MRI  protocols  capture  different  aspects  of  the  Ms  disease  puzzle.  We  have  therefore  taken  a  multiprotocol  approac 
herein  we  uEaS  major  protocols  currently  employed  in  MS  imaging  and  attempt  to  get  a  °f 

th*  most  snedfic  wav  of  characterizing  the  disease  through  image  derived  parameters.  In  Section  we  .  . 

™oCd1toStong tomSL  from  toe  various  MRI  protocols.  and  in  Section  3.  we  sommarize  the  toicaltori 
“EtSSTSe  to  sS  4.  we  summarize  to,  lessons  toed  from  tois  ezperienc.  ami  toe  omamndmg  problems 

that  need  to  be  addressed. 

2.  METOHDS 

2.1  Data  acquisition 

The  MRI  protocols  we  have  been  using  in  our  patient  studies  are  listed  in  the  following  table.  Every  patient  recruited  in  our 
S^^aSS?on7Sted  in®  the  table.  Our  image  database  currently  consists  of  die  fofiowmg jnumte -of  pg* 
studies  and  3D  volume  images  that  have  been  processed  by  the  methods  described  m  660 

studies  each  with  and  without  contrast,  670  MT  studies,  altogether  100  patients  and  over  4000  3D  scenes. 


Pulse  Seq. 


FSEVE 


Plane 


AXIAL 


3DMT  Vas  Tof  I  AXIAL 


SE _ 

SE  with  Gd 


AXIAL 

AXIAL 


TR  (ms) 


2500 


106 


600 

600 


TE  (ms) 


18,90 


5 


27 

27 


Slice  Th  (mm) 

Matrix 

NEX 

FOV  (cm) 

3 

192 

1 

22 

5 

128 

1 

22 

3 

192 

1 

22 

3 

192 

1 

22 

The  T2  and  PD  images  in  these  acquisitions  are  in  registration  since  they  are  acquired  simultaneously.  Sin“^^ 

images  are  also  in  registration.  However,  between  these  two  sets  and  among  other  image  data,  registration  cannot  be 

guaranteed. 

2JL  Fuzzy  connectedness  image  segmentation 

The  method  of  fuzzy  connectedness  forms  the  essential  underpinning  of  all  segmentation  algorithms  utilized  m  our  MS 
image  analysis  efforts.  Therefore,  we  will  first  give  an  outline  of  its  principles.  These  principles  are  applicable  to  n- 
dimensionaT vector- valued  scenes.  However,  our  description  will  be  confined  to  the  three-dimensional  case  and  to  scalar 

valued  scenes. 

We  represent  a  volume  image,  called  a  3D  scene  (a scene  for  short)  6  by  a  pair  e  =  (C.  f)  where  C  is  a  rectangular  array  of 
voxels,  called  the  scene  domain ,  and/is  a  function  that  assigns  to  every  voxel  ceCan  intensity'  value^from  a  raige  [LH] 
where  L  and  H  are  integers.  Objects  such  as  WM  in  the  brain  are  manifest  in  scenes  with  a  heterogeneity  of  property  values 


(such  as  PD)  because  of  the  heterogeneity  of  material  composition  inherent  in  the  object,  and  noise,  blurring  and  background 
variations  introduced  by  the  imaging  device.  In  spite  of  this  graded  composition  of  intensity  values  within  object  regions, 
human  readers  perceive  regions  in  the  scene  belonging  to  the  same  object  as  a  whole  (gestalt).  This  property  of  hanging- 
togetherness  of  image  elements  and  their  graded  composition  are  both  fuzzy  phenomena,  which  should  be  addressed 
properly  for  effective  image  segmentation.  While  the  property  of  graded  composition  has  been  handled  in  the  past  through 
fuzzy  logic  and/or  probabilistic  methods,  fuzzy  connectedness  was  the  first  framework  that  allowed  capturing  the  idea  of 
fuzzy  hanging-togetherness  via  fuzzy  topological  principles.32  It  is  defined  as  follows. 

We  think  of  any  nearby  voxels  c  and  d  in  C  as  having  a  fuzzy  adjacency.  The  strength  )Jdc,  d)  of  this  fuzzy  relation  a  is 
smaller  the  farther  c  and  d  are.  The  idea  behind  a  is  to  capture  the  blurring  property  of  imaging  devices.  We  define  another 
fuzzy  relation  Kin  C  that  assigns  to  every  pair  of  voxels  (c,  d)  an  affinity.  The  strength  ji^c,  d)  of  this  relation  depends  on 
how  near  c  and  d  are  spatially  (i.e.,  on  pjc,  dj)  and  on  how  similar  the  scene  intensities  f[c)  dndfid)  are  as  well  as  how 
similar  intensity-based  properties  computed  at  c  and  d  are.  The  intent  here  is  that  k  is  "local”;  that  is,  if  c  and  i  are  far  apart 
(spatially),  then  their  affinity  is  0.  To  define  fuzzy  connectedness  K  as  another  fuzzy  relation  in  C,  we  consider  all  possible 
"paths"  between  all  possible  pairs  of  voxels  in  C.  A  path  between  any  voxels  c  and  d  is  simply  a  sequence  of  nearby  voxels 
starting  from  c  and  ending  in  d.  To  every  path  from  c  to  d  we  assign  a  "strength  of  connectedness"  which  is  simply  the 
smallest  affinity  of  pairwise  voxels  along  the  path  (weakest  link).  Finally,  the  strength  of  connectedness  }J^[ct  d)  between  c 
and  d  is  the  largest  of  the  strengths  of  connectedness  of  all  paths  between  c  and  d.  In  determining  a  fuzzy  connected  object  in 
C,  jJLidc,  d)  should  be  determined  for  all  possible  pairs  (c,  d)  of  voxels  in  C.  This  computationally  explosive  problem  is 
considerably  simplified  through  some  key  theoretical  results,  and  finding  fuzzy  connected  objects  essentially  reduces  to 
dynamic  programming.32 

The  definition  of  affinity  is  a  key  to  the  effectiveness  of  segmentation  in  this  method.  We  think  of  the  strength  of  affinity 
pUc,  d)  between  c  and  d  to  consist  of  three  components:33 

uJc-d) =*(w"*)  0> 

Ha  is  the  adjacency  component  mentioned  above.  nr  is  a  homogeneity-based  component.  As  per  this  component,  the  greater 
the  homogeneity  of  the  intensity  in  the  vicinity  of  c  and  d  is  the  greater  is  the  affinity.  tu  is  an  object-feature-based 
component.  As  per  this  component,  the  closer  the  intensity-based  features  at  c  and  d  are  to  an  expected  value  of  these  features 
for  the  object  region,  the  greater  is  their  affinity.  A  functional  form  of  g  we  have  used  commonly  is 

Hk  ( c,d)  =  va  (c,  d (c,  d)  (c,  d)  .  (2) 

In  the  scale-based  approach,33,5*  d)  and  /Vc»  d)  are  defined  taking  into  consideration  all  voxels  within  a  sphere  centered 
at  c  and  d.  The  radius  of  this  sphere  is  related  to  the  "scale"  at  c  and  d  of  the  object  being  defined.  The  scale  of  the  object  at 
any  voxel  c  is  defined  as  the  radius  of  the  largest  sphere  that  can  be  placed  with  its  center  at  c  such  that  it  encloses  only 
voxels  in  the  object  region.  Paradoxically,  this  may  sound  like  we  need  object  segmentation  to  define  scale  at  c.  We  have 
developed  a  simple  algorithm31 33  that  does  not  require  object  segmentation  for  estimating  scale  at  every  voxel  c.  The 
algorithm  uses  simply  a  homogeneity  measure  to  determine  at  what  radius  there  is  a  sudden  change  in  homogeneity  as  the 
sphere  is  increased  from  a  radius  of  1.  We  have  shown33  that  scale-based  fuzzy  connectedness  is  less  sensitive  to  noise  and 
can  detect  thin  and  subtle  aspects  of  the  object  more  effectively  than  the  original  fuzzy  connectedness  method,  yet  the  theory 
of  the  original  fuzzy  connectedness  framework  still  holds.  Scale-based  fuzzy  connectedness  is  a  powerful  segmentation 
method  that  differs  from  published  methods  in  that,  (1)  it  takes  into  account  in  its  design  the  size  of  the  object  in  different 
parts  of  the  scene  and  image  noise  and  other  artifacts;  (2)  it  handles  both  the  graded  composition  and  hanging-togetherness  in 
a  fuzzy  setting.  We  note  that  since  the  local  scale  is  taken  into  account  in  determining  d),  the  actual  functional  form  of 
the  affinity  varies  over  the  scene  domain,  adapting  to  the  local  object  size  and  noise  characteristics.  We  have  also  shown  that 
slow  background  variations  do  not  affect  fuzzy  connected  object  segmentations.32,33 

In  the  rest  of  this  section,  we  summarize  our  approach  to  analyzing  the  scenes  for  the  different  protocols.  The  basic  premise 
behind  these  approaches  is  the  following.  We  think  of  segmentation  to  consist  of  two  related  tasks  -  recognition  and 
delineation.  Recognition  is  the  high-level  process  of  roughly  determining  the  whereabouts  of  the  object  in  the  scene. 


n  iin*atu>n  is  the  low-level  process  of  determining  the  precise  spatial  exrat  and  voxel-by- voxel  material  percent  content  of 
J^JL,  Knowledgeable  humans  can  outperform  computer  algorithms  in  most  recognition  tasks  such  as  diose  encountered 
*eSSn Smely,  computer  algorithms  can  be  devised  to  more  precisely,  accurately.  and  efficiently  dehneme 
!i2s?scenes  than  manual  delineation.  Clearly,  manual  delineation  to  indicate  the  graded  material  composition  wuhman 
"KTis  Sssible.  The  system  we  have  developed  for  the  above  four  tasks  exploits  this  synergy  between  human  and 
Snputer  abilities  in  devising  practical  solutions  that  can  be  used  routinely  in  clinical  trials. 


2J.  Analysis  of  FSE  T2,  PD  images 


The  approach  here  consists  of  the  following  four  steps.33  We  consider  here  the  scene  to  be  vector-valued  with  two  (T2  and 
pD)  values  per  voxel. 


1. 


On  one  slice,  roughly  centrally  situated  in  the  brain,  an  operator  specifies  a  few  voxels  (seeds).  This  is  a 
recognition  step  that  uses  the  superior  human  knowledge.  Points  (voxels)  are  specified  for  the  CSF,  G  ,  an 
WM  regions  (and  not  for  lesions). 


2.  This  is  a  delineation  step.  The  seed  voxels  are  utilized  to  determine  the  fuzzy  connected  objects  that  contain 
This  results  first  in  a  segmentation  of  WM,  GM,  and  CSF.  This  knowledge  is  subsequently  used  to 
determine  automatically  a  set  of  points  (determined  as  holes  in  the  union  of  GM  and  WM  fuzzy  objects)  in 
each  3D  object  that  are  used  subsequently  to  delineate  the  lesions,  each  as  a  3D  fuzzy  connected  object 


3. 


This  is  again  a  recognition  step  taking  help  from  a  human  operator.  In  this  step,  the  operator  accepts  true 
Itams  and  rejects  false  lesions  with  a  mouse  button  click.  Each  3D  lesion  is  displayed  on  one  slice  image 
passing  close  to  the  centroid  of  the  lesion.  The  operator  may  override  this  mode  of  display  and  examine  the 
ipcinn»  on  all  or  any  selected  slices.  All  false  positives  are  eliminated  in  this  fcshion.  These  are  usually 
arti&cts  and  choroid  plexus.  We  found  that  the  number  of  false  negatives  is  very  low  in  our  system  (see  next 
section  for  H»»ii<  on  validation).  Nevertheless,  the  system  allows  in  this  step  selecting  new  seed  points  for 
lesions.  These  are  subsequently  utilized  in  delineating  these  lesions  by  repeating  part  of  Step  2  above. 


4  The  final  step  is  the  computation  of  quantitative  measures  from  the  segmented  objects.  These  includethe 
number  of  3D  lesions  and  their  total  volume  (T2LV),  the  volume  of  CSF  (CSFV),  the  volume  of  the  brain 
parenchyma  (BPV)  which  is  the  volume  of  the  union  of  GM  and  WM,  and  normalized  BPV,  nB 
BPV/(BPV  +  CSFV).  The  purpose  of  the  last  measure  is  to  express  the  parenchymal  volume  independent  of 
subject-to-subject  variations  in  brain  size. 


2.4  Analysis  of  TIWIE  images 

The  approach  here  consists  of  the  following  steps.34 

1.  A  conservative  threshold  is  determined  automatically  from  the  histogram  of  the  TIWIE  3D  scene  (see  for 
details).  The  purpose  of  this  threshold  is  to  select  seed  points  automatically  within  die  enhancing  lesions. 
This  is  a  recognition  step. 

2.  The  seed  points  are  utilized  to  determine  the  fuzzy  connected  objects  that  contain  them.  This  delineation  step 
often  results  in  the  delineation  of  vessels.  Taking  the  volume  enclosed  by  the  object  as  a  criterion,  large 
objects  (large  vessels)  are  automatically  discarded. 

3.  This  step  is  identical  to  Step  3  of  the  previous  section. 

4.  For  the  accepted  lesions,  their  total  volume  (T1EV)  and  number  (nTlEV)  are  computed. 


2.5.  Analysis  of  MT  images 

The  steps  involved  in  this  process  are  as  follows.37 


I.  The  torn  is  firs.  segmented  tan  the  two  MT  scttnes  fcorresponding  «,  dte  off-resonance  puise  on  end  off, 

The  train  mask  is  then  utilized  to  comoute  an  MTR  scene  a"  -( C *  f*\  j,-  ,  . 

«  ‘  -lc«7 *-4J.  which  1S  such  that 

^  ^  M  will  have  a  vaiue  as  detomined 

from  the  two  MT  scenes  as  described  in  Section  1.3. 

1  ZSXZSSZSXSSS TZ:£2?G£X=! 

segmented  for  WM  ami  GM  as  described  in  Section  13.  Tie  MTR  scenes  ,«  and 

subsequently.  S  are  computed 


preach  of  the  ,vrre  scenes  <?'  eZ .  and  .  a  btstogmm  of  die  whole  scene  normalised  to  BPV  it 
WM,  or  GM.  peaCwntiles  (P-5,,  P50*)  and  the  mean  m„  where  x  stands  for  BP, 


3.  RESULTS 


3.L  Validation 


effort  to  mate  the 

A'cmacy  denotes  the  degri  of  agreeS* dSSrtS  Wee**  input  required  by  the  algmitos. 

(the  degree  of  automation),  which  may  be  expressed  as  the  ooeratOT  1"cJcaIes  *e  degree  of  operator  help  required 

is  usually  required  to  achieve  an  ImprLnen^oSf B A  C°mpr0mise  111  « 
vanous  analysis  processes  described  in  the  previous  section  We  fc^Td^t?  .prCCIS1°n'  accuracy>  and  efficiency  of  the 
account  the  subjectivity  that  appears  in  Step?  1  and  3  and  in  the  Dlacim^r  nf  ^  evalua“ons  of  our  system  taking  into 
results  is  as  follows.35"37, 39  For  T2JLV*  Infer  jmH  inm  ^  ^atLcnt  m  scanner.  A  summary  of  the 

1.5*1  95%  confidence  intiS  SgSve^TS?  °f  Var“°"  '  °'9*  "f“  »•  - 

operator  variations,  and  1.3%  false  negative  TlELvf^ean  ^  **  SparC'  F°r  T1ELV:  No  inter-  and  intra- 

operator  coefficient  of  variation  <  0.5%includinz  TL  **  StUdy'  For  BPV:  and  inn- 

£=££  MT* » "• 

—  -  to  .ha. 

tolttae  Bkmg  tecognidcn  help  tan ‘the  user  to  nSSTfiZ  ST?"  «?“  “•  by  design.  not  autcanadc  tn 
efficaciously  as  possible  to  make  the  sysmTSfi^d H°w'ver’  ““  “P  *  taken  as 
without  sacrificing  precision  and  accuracy,  ^  '  e  continue  to  improve  the  efficiency  of  our  system 

3JL  Clinical 


clinical  conelaticn^J^rith  dK^sTscor^0^^  “  ^  “  analy7ing  the  ima§es  h  our  database  and  in  conducting 

group  and  ^ 0^)^-  ^  "BPV  ('0-73’  ^  <  -001)  in  the  RR 

wtth  nBPV  (-0  66,  p  <  0.02)  £  the  CP^oup M  ^  SSon  ^  a  good  c“n 

Generally  we  found  CSFV  to  be  significantly  larger  in  MS  patients  and  nRPV  J*1*  300123  was  not  ^gnificant. 

control  subjects  ( p  <  0.005).  A  good  correlatio,  seen  IwKvi^ffl  ^  *  age'matched  nonnal 

lesions, .«  found  a  sipifiesn. »  <  0.05,  difference  m  toe  of  *e  U  nemopsychologii' «  S£SS 


and  without  these  lesion is,  possibly  indicating  that  %of*e  toal  T2LV  while  deep  gray 

lesions  in  gray  maner,  we  found  drat  ccudcal  Spy  ^  found  between  these  volumes  and  clinical  disabiiiQt.  By 
Inane, lemons  compnsedahootu.6%.  No  agrafe* ««"*«“ ^^"^eladcn  widi  disability  was  nBPV  (-0.69 ,p< 

^  -V  be  an  imponant  measrae  to  characranrang  draease  * 

X^WI. 

significant  decrease  in  the  number  of  enhancing  lesions  (p  -  •  )  dssue  ^  significantly  less  {p  =  0.02)  in  the  group 

No  differences  were  seen  ®J2LV  in  the  wo  gr°“PS^  w  report  die  MRI  activity  in  copaxone  clinical  trials  of 

^jing  copaamte  tot  ^ugh  ^  MS  software  sysram  dtat  the  drag  made  a  significant  unpac.  on 

inflammation  and  on  preventing  brain  atrophy  in  RR  patients. 

MT  imaging:^3  4*-50  Compared  to  the  normal  control  nwS 

patients)  was  significandy  lower  ip  <  0.0005),  and  so  were  Map  ip  -  ^  parameter  correlated 

Ocularly,  showed  the  best  correlation  ^  witi  CSFV^S, p <  0.0001) 

Songly  with  nBPV  in  both  RR  (0.69,  p  <  0.001)  and  theCP  (0.87,  p  <  auui;  ^  duration  ^  <  0.0l). 

when  the  two  groups  are  pooled  into  a  angle  gro“P-  (_  <  q  001).  The  unnocmalized  histogram  peak  height 

Individual  neuropsychological  teas  correlate  wnhMTR  measures  ip  <  0^.  ^  a  serial  MW  study  of  the  RR 

differed  ip  <  0.05)  among  severely  impaired,  moderamly  ^  ^ile  no  significant  changes  in  EDSS  were 

group.  nPHa,  demonstrated  a  subde  but  ag^cau £  < *an  in  nSmal  controls."  In  longitudinal 
noted.  Mwm  and  mwM  significantly  .  5)  increasing  diCM<a*  duration.  Up  to  44%  of  new  lesions 

studies,  both  these  entities  shifted  in  the  Section  o ^  M  abn0rnial  by 

identified  on  later  studies  were  demonstrated  to  J  (D<0  01)  lower  in  RR  MS  patients  than  in  normal  controls. 

MIR  criteria.  As  to  the  GM,  Mgm  and  mo*  were  sigrufic^ntly  (p  <  0.01)  lower  in  **  y 

nPHcsM  inversely  correlated  with  EDSS  (r  =  0.65,  p  <  0.01). 

Insorararay.  ®c  surges,— 

d  CONCLUDING  REMARKS 


4.L  Automation,  failure,  user  assistance 

My  a*raenadra,  racdrad  cun,  . d  M.  will.  go  -hJB^^ 

clinical  Ural  mding.  Tteetoe.  to  qualuy  hhtounce.  »  »  ‘p^  rJdrabic  after  only  »  great  deal  of 

processing  loop.  Complete  automation  may,  therefor  ,  ^  researchers  and  developers,  our  aim 

experience  within  the  same  imaging  modality,  Pro*“*?  ’  .  ^  segmentation  algorithms  with  a  research 

should  be  to  consider  human  interaction  within  *e  2 <LShere  may  not  be  optimal  yet  in  terns 

goal  of  minimizing  this  interaction  as  much  as  possible.  Die  system  dia^ mw  ^  efficiency, 

of  the  degree  of  automation  achievable,  but  it  is  certainly  practical.  We  continue  to  tighten 

4.2.  Standardization  of  MR  intensity  scale 

A  major  difficulty  with  the  MRI  techniques  for  most .Prowcols  has^ *Tsln^l£Sl 
even  within  the  same  protocol,  for  the  same  body  region,  or  g  importantly,  in  image  segmentation  and 

poses  problems  in  finding  the  proper  window  setting  for  J ? protocol-specific  intensity  meaning,  setting 

as— 


We  have  recently  developed  a  method*3,  u  to  standardize  the  MR  image  intensity  scale.  It  is  a  post-acquisition  processing 
method  which  maps  nai-Iineariy  the  scene  intensities  in  any  given  scene  acquired  as  per  a  given  protocol  for  a  given  body 
region  into  a  standard  scale  so  that  the  same  intensities  in  the  transformed  scene  will  have  the  same  tissue  specific  meaning. 
The  method  is  based  on  deforming  histograms  so  that  they  are  as  similar  as  possible  for  all  scenes  of  the  same  body  region 
and  protocol.  We  are  in  the  process  of  utilizing  this  method  in  modifying  our  system  to  improve  the  efficiency  of  the  various 
processes  of  quantification  without  sacrificing  precision  and  accuracy. 

4J.  Standard  intensity-based  analysis 

All  objects  -  normal  and  pathological  -  have  a  heterogeneous  tissue  composition.  This  combined  with  the  blurring  and  noise 
introduced  by  the  imaging  device  makes  tissue  regions  have  a  heterogeneity  of  intensities.  We  believe  that  this  heterogeneity 
has  tissue-specific  information  and  is  useful  in  characterizing  disease  stage  and  severity.  Such  a  characterization  becomes 
impossible  when  there  is  scanner-dependent  intensity  variation.  As  illustrated  by  MTR  analysis,  because  of  the 

tissue-specific  meaning  of  MTR  values,  standardization  may  permit  us  to  treat  other  protocols  (T2,  PD,  T1WIE,  T1WI)  also 
in  the  same  tissue-specific  way  as  MT  imaging.  More  importantly,  flat  measures  such  as  volume  of  tissue  regions  ignore  the 
heterogeneity  information  and  may  lose  important  disease  specific  information.  We  should  really  consider  volume 
distributions,  that  is,  intensity  histograms  within  segmented  tissue  regions,  as  demonstrated  by  MTR  analysis,  for 
understanding  subtle  disease  processes. 

4.4.  MS  segmentation  "workshop" 

hi  MS  (and  ocher  neurological  applications),  a  variety  of  MRI  protocols  are  utilized  as  we  examined  in  this  paper.  The  actual 
imaging  parameters  used  in  these  protocols  vary  among  institutions.  In  spite  of  the  numerous  brain  MR  image  segmentation 
methods  developed  during  the  past  15  years,  none  of  them  is  capable  of  handling  variations  within  the  same  protocol,  and 
much  less,  the  variations  among  protocols.  What  we  need  is  a  segmentation  "workshop"  wherein  a  protocol-specific 
segmentation  method  can  be  quickly  fabricated.  For  the  MS  application,  we  believe  that  the  fuzzy  connectedness  framework 
can  be  utilised  to  build  such  a  workshop  and  we  are  working  toward  this  goal. 
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ABSTRACT 

In  the  past,  we  have  presented  three  user-steered  image  segmentation  paradigms:  live  wire,  live  lane,  and  the  3D 
extension  of  the  live-wire  method.  In  this  paper,  we  introduce  an  ultra-fast  live-wire  method,  referred  to  as  iive- 
wire-on-the-fly,  for  further  reducing  user’s  time  compared  to  live  wire.  For  both  approaches,  given  a  slice  and  a 
2D  boundary  of  interest  in  this  slice,  we  translate  the  problem  of  finding  the  best  boundary  segment  between  any 
two  points  specified  by  the  user  on  this  boundary  to  the  problem  of  finding  the  minimum-cost  path  between  two 
vertices  in  a  weighted  and  directed  graph.  The  entire  2D  boundary  is  identified  as  a  set  of  consecutive  boundary 
segments,  each  specified  and  detected  in  this  fashion.  A  drawback  in  live  wire  is  that  the  speed  for  optimal  path 
computation  depends  on  image  size,  compromising  the  overall  segmentation  efficiency.  In  this  work,  we  solve  this 
problem  by  exploiting  some  properties  of  graph  theory  to  avoid  unnecessary  minimum-cost  path  computation  during 
segmentation.  Based  on  164  segmentation  experiments  from  an  actual  medical  application,  we  demonstrate  that,  live- 
wire-on-the-fly  is  about  1.5  to  33  times  faster  than  live  wire  for  actual  segmentation,  although  the  pure  computational 
part  alone  is  found  to  be  over  a  hundred  timx  fester. 

Keywords:  image  segmentation,  boundary  detection,  active  boundaries,  3D  imaging,  shortest-path  algorithms, 
dynamic  programming,  graph  theory. 


1.  INTRODUCTION 

Image  segmentation  is  a  hard  problem  with  numerous  applications  in  the  imaging  sciences.1  It  consists  of  two  tightly 
coupled  tasks  -  recognition  and  delineation.  Recognition  is  the  process  of  identifying  roughly  the  whereabouts  of  a 
particular  object  in  the  image  and  delineation  is  the  process  of  specifying  the  precise  spatial  extent  and  composition 
of  this  object.  While  computer  algorithms  are  very  effective  in  object  delineation,  the  absence  of  relevant  global 
object-related  knowledge  is  the  main  reason  for  their  failure  in  object  recognition.  On  the  other  hand,  a  simple  user 
assistance  in  object  recognition  is  often  sufficient  to  complement  this  deficiency  and  to  complete  the  segmentation 
process.  There  are  many  difficult  segmentation  tasks  that  require  a  detailed  user  assistance.  To  address  these 
problems,  a  variety  of  interactive  segmentation  methods  are  being  developed.2  These  methods  range  from  totally 
manual  painting  of  object  regions  or  drawing  of  object  boundaries  to  the  detection  of  object  region/boundaries  with 
minimal  user  assistance.3^7 

We  have  been  developing  interactive  segmentation  strategies  with  two  specific  aims:  (i)  to  provide  as  complete  a 
controi  as  possible  to  the  user  on  the  segmentation  process  while  it  is  being  executed,  and  (ii)  to  mini™;,,,  the  user 
involvement  and  the  total  user’s  time  required  for  segmentation,  without  compromising  the  precision  and  accuracy  of 
segmentation.  Our  strategy  in  these  methods  has  been  to  actively  exploit  the  superior  abilities  of  human  operators 
(compared  to  computer  algorithms)  in  object  recognition  and  the  superior  abilities  of  computer  algorithms  (compared 
to  human  operators)  in  object  delineation. 

In  the  past,  we  have  presented  two  user-steered  segmentation  paradigms,  referred  to  as  live  wire  and  live  lane,6-8  to 
segment  3D/4D  object  boundaries  in  a  slice-by-slice  fashion.  These  methods  are  in  routine  use  in  two  applications8-12 
with  over  15,000  tracings  done  so  far.  Although  the  live- wire  method  has  its  origin  in  some  early  joint  work  between 
Barrett  and  Udupa, 13-15  this  method  has  been  subsequently  developed  independently  by  the  two  groups.6-8,18-18 
There  are  many  differences  between  the  live-wire  method  developed  by  each  group,  as  previously  explained  in.6 
Besides  these  differences,  we  have  extended  the  ideas  underlying  the  live-wire  method  to  create  new  methods,  live 
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lane6  and  the  3D  extension  of  live  wire.16  In  this  paper,  we  introduce  an  ultra-fast  live-wire  method,  referred  to 
as  live-wire-on-the-fly,  with  a  new  live-wire  algorithm  for  drastically  reducing  user’s  time  compared  to  our  previous 
work  on  2D  live  wire.  F 

In  live  wire,6-7  to  segment  a  2D  boundary,  the  user  initially  picks  a  point  on  the  boundary  and  all  possible 
minimum-cost  paths  from  this  point  to  all  other  points  in  the  image  are  computed  via  dynamic  programming. 
Subsequently,  a  “live  wire”  is  displayed  in  real  time  from  the  initial  point  to  any  subsequent  position  t.alron  by  the 
cursor.  If  the  cursor  is  close  to  the  desired  boundary,  the  live  wire  snaps  on  to  the  boundary.  The  cursor  is  then 
deposited  and  a  new  live-wire  segment  is  found  next.  The  entire  2D  boundary  is  specified  via  a  set  of  live-wire 
segments  m  this  fashion.  A  drawback  of  this  approach  is  the  computational  time  for  all  possible  minimum-cost 
segments  from  each  selected  point  on  the  boundary  to  other  points  in  the  image.  This  time  increases  with  the  size 
of  the  image  compromising  the  interactivity  of  the  method  in  some  practical  situations.  For  images  from  256  x  256 
to  1024  x  1024  pixels,  for  example,  live  wire  running  on  a  300MHz  Pentium  PC  requires  about  2  to  180  seconds" to 
compute  all  possible  minimum-cost  segments  from  each  selected  point. 

In  live  wire  on  the  fly,  the  user-interaction  process  remains  the  same,  but  we  have  devised  a  linear  time  complexity 
algorithm  to  save  a  considerable  amount  of  user  time  by  avoiding  the  computation  of  all  possible  minimum-cost 
segments.  When  the  user  selects  a  point  on  the  boundary,  the  live-wire  segment  is  computed  and  displayed  in  real 
time  from  the  selected  point  to  any  subsequent  position  of  the  cursor  in  the  image.  To  make  this  feasible,  we  exploit 
the  fact  that  by  the  time  we  have  found  a  live-wire  segment  with  cost  value  K,  we  have  actually  found  all  possible 
live-wire  segments  with  cost  value  less  than  K  in  the  image.  Moreover,  any  live-wire  segment  with  cost  value  greater 
than  or  equal  to  K  contains  one  of  the  previous  live-wire  segments  with  cost  value  less  than  K.  Therefore,  the 
computation  of  the  minimum-cost  segment  from  a  selected  point  to  the  current  position  of  the  cursor  uses  the  results 
of  computation  from  the  selected  point  to  the  previous  position  of  the  cursor. 

In  Section  2,  we  present  the  live-wire-on-the-fly  method  and  its  algorithms.  In  Section  3,  we  present  the  results 
of  evaluation  between  live  wire  and  live  wire  on  the  fly  based  on  efficiency  for  segmentation.  Finally,  we  state  some 
concluding  remarks  in  Section  4. 


2.  LIVE-WIRE-ON-THE-FLY 

We  define  a  2D  scene  C  as  a  pair  (C,g)  consisting  of  a  finite  2D  rectangular  array  C  of  pixels  and  a  function 
:  C  -4  [L,  H]  that  assigns  to  each  pixel  p  in  C  an  intensity  value  lying  in  an  interval  [I,  H].  We  associate  with 
C  a  directed  graph  m  which  the  vertices  of  the  pixels  represent  the  nodes  of  the  graph  and  the  oriented  pixel  edges 
represent  the  arcs.  A  2D  boundary  of  interest  in  C  is  a  closed,  oriented,  and  connected  contour  made  up  of  oriented 
pixel  edges.  Each  oriented  pixel  edge  in  C  is  a  potential  boundary  element  6.  which  is  called  a  bel  for  short  To  each 
bel  6,  we  assign  a  set  of  features  whose  values  characterize  the  “boundariness”  of  b.  These  values  are  converted  to 
a  single  joint  cost  value  c(b)  per  bel  b.  Thus,  the  problem  of  finding  the  best  boundary  segment  (live-wire  segment) 
between  any  two  points  (pixel  vertices)  specified  on  the  boundary  is  translated  to  finding  the  minimum-cost  path 
between  the  corresponding  two  vertices  of  the  graph.  The  issues  about  selection  of  features  and  how  to  convert 
feature  values  into  cost  values  were  previously  addressed  in.6  The  problem  we  want  to  address  here  is  how  to  reduce 
tne  tune  for  optimum  path  computation,  and,  consequently,  the  total  user’s  time  required  for  segmentation. 

To  tackle  this  problem,  we  will  exploit  some  known  properties  of  graph  theory,  particularly  for  the  computation 
of  shortest-paths,  as  described  in  Section  2.1.  This  leads  to  the  algorithms  presented  in  Section  2.2. 

2.1.  Graph  Properties  of  Shortest  Paths 

In  the  literature  on  shortest-path  algorithms,19  there  are  many  efficient  solutions  for  finding  minimum-cost  paths  in 
wag  ed  and  directed  graph.  Particularly,  we  have  adopted  Dial’s  implementation  of  the  Dijkstra’s  algorithm  20 
1  his  algorithm  computes  the  shortest-paths  to  all  nodes  from  a  single  node  in  0{m  +  nC)  time,  where  m  is  the 
number  of  arcs,  n  is  the  number  of  nodes,  and  C  is  the  maximum  cost  assigned  to  any  arc  in  the  graph.  Actually 
“  this  case  the  cost  assigned  to  each  arc  should  be  an  integer  in  the  interval  [0.CJ.  Dial’s  solution  uses  a  circuit 
queue  with  C  +  1  buckets  of  nodes  as  the  priority  queue  of  the  Dijkstra’s  algorithm.  Since  the  bottleneck  of  the 
Dijkstras  algorithm  is  in  maintaining  the  priority  queue,  Dial’s  solution  uses  the  bucket  sort  algorithm  to  speed  up 
this  process.  We  will  come  back  to  this  issue  in  Section  2.2.  y  H 

In  our  problem,  the  live-wire  segment  between  a  selected  point  v,  on  the  boundary  and  the  current  position 
of  the  cursor  in  C  is  the  shortest-path  P  =  (v9  ^  ve)  from  v,  to  vc  in  our  graph,  where  the  cost  of  P,  denoted 


K(P)y  is  the  sum  of  the  joint  costs  c(b)  of  all  bels  b  comprising  P.  In  fact,  Dijkstra’s  algorithm  returns  a  tree  of 
minimum-cost  path  (or  a  tree  of  shortest-path)  rooted  at  v$  ,21  which  consists  of  all  minimum-cost  paths  from  vt  to 
all  vertices  in  C.  We  will  denote  this  tree  by  T(vs)  =  {P  =  ( va  ve)/v€  €  C}. 

For  any  real  number  fc,  we  denote  by  Tk(vs)  the  tree  of  minimum-cost  path  rooted  at  vs  such  that  the  cost  of 
any  path  in  this  tree  is  less  than  Jfc.  That  is,  Tk{v ,)  *  {P  =  {vs  —  ve)/ve  €  C,K{P)  <  k}.  The  algorithm  reported 
in  this  paper  exploits  the  following  properties  of  T(v,). 

1.  To  compute  the  minimum-cost  path  P  =  ( vs  ve)  with  cost  K(P ),  there  is  no  need  to  compute  Tk(vs)  for 

k  >  K{P). 

2.  By  the  time  we  have  found  the  minimum-cost  path  P  =  {vs  —  t/e)  with  cost  K{P ),  we  have  actually  found  the 
tree  of  minimum-cost  path  Tk(p){vs). 

3.  The  tree  of  minimum-cost  path  Tk{vs)  contains  the  tree  of  minimum-cost  path  TK(p){vs)  whenever  k  >  K{P). 

We  use  the  first  property  to  modify  Dial’s  implementation  of  the  Dijkstra’s  algorithm  to  quit  optimum  path 
computation  by  the  time  we  have  found  the  minimum-cost  path  P  =  (vs  ^  ve).  We  call  this  algorithm  DSP  (see 
Section  2.2).  We  use  the  second  property  to  avoid  optimum  path  computation  for  any  path  P '  =  (nj  ***  vfe)  with 
cost  K(P‘)  <  K(P).  Thus,  when  the  user  moves  the  cursor  to  a  new  position  v'e ,  such  that  K(P')  <  K{P ),  and  we 
have  already  found  P,  the  algorithm  just  shows  P'  =  ( v9  ^  v'c )  without  requiring  computation.  We  use  the  third 
property  to  continue  optimum  path  computation  for  paths  Pf  =  ( vs  vfe)  with  costs  K{Pt)  >  K(P)  based  on  the 
previous  result  of  algorithm  DSP  for  computing  P. 

2.2.  ALGORITHMS 
Algorithm  LWOF 

Input:  The  joint  cost  function  c  and  an  initial  vertex  vo  selected  on  a  2D  boundary  of  interest  in  C; 

Output:  A  closed,  connected,  and  oriented  contour  B  (made  up  of  bels); 

Auxiliary  Data  Structures:  A  2D  “cumulative  cost”  array  cc  representing  the  total  cost  of  the  optimal  paths  found 
so  far  from  vs  to  other  vertices  in  C;  a  2D  “direction”  array  dir  indicating,  for  each  vertex,  to  which  of  its  immediate 
neighboring  vertices  the  optimal  path  goes;  a  circular  queue  Q  of  vertices  with  C  4*  1  buckets;  a  list  L  of  vertices  which 
have  already  been  processed;  a  current  path  P(vs  ve),  where  vs  is  the  current  point  selected  on  the  boundary  and 
ve  is  the  current  position  of  the  cursor  in  C;  and  a  list  B  of  bels  which  have  already  been  identified  as  belonging  to 
the  boundary  of  interest  in  C; 


begin 

1.  set  cc(v)  to  oo  and  dtr(t;)  to  null  for  all  vertices  v  in  C,  and  set  L  to  empty; 

2.  vs  <-  v0 ,  set  cc(vs)  to  0,  and  put  vs  in  Q; 

3.  repeat 

a.  determine  the  vertex  ve  in  C  pointed  to  by  the  cursor; 

b.  if  vt  is  not  a  vertex  of  any  bel  in  B  then 

(i)  compute  P  i-DSP(u5, uc,Q,cc,c, dir, L)  and  display  the  bels  in  P; 

(ii)  if  vt  is  selected  by  the  user  and  vt  €  C  then 

a.  add  the  bels  in  P  to  B; 

b.  set  cc(u)  to  oo  and  dir(v)  to  null  for  all  vertices  v  in  C; 

c.  remove  all  vertices  v  from  Q,  and  remove  from  L  all  vertices  v  which  do  not  belong  to  any  bel  in 

B; 

d.  v3  <-  ve,  set  cc(vs)  to  0,  remove  vs  from  L,  and  put  u,  in  Q\ 


endif 

endif 

until  the  user  indicates  a  “closer  operation; 

4.  ve  4-  vo  and  remove  vt  from  I; 

5.  compute  P  *-DSP (vSyvt,Q,cc,c,dir,L)  and  display  the  bels  in  P\ 

6.  add  the  bels  in  P  to  B  and  output  the  bels  in  B; 

end 


Algorithm  DSP 


Input:  an  initial  vertex  v9\  a  terminal  vertex  vt\  the  circular  queue  Q\  the  cumulative  cost  array  cc;  the  joint 
cost  function  c;  the  direction  array  dir ;  and  the  list  L  of  already  processed  vertices; 

Output:  A  set  of  bels  forming  an  optimal  path  from  vt  to  vt\ 


begin 

1.  while  ve  £  L  do 

a.  remove  a  vertex  v  from  Q  such  that  cc(u)  =  miiVeQ{cc(t/)},  and  put  t;  in  L; 

b.  for  each  vertex  v*  such  that  v*  is  in  the  set  of  the  4-adjacent  neighbors  of  v  and  v'  £  L  do 

(i)  compute  cctmp  =  cc{v)  -f  c(V)  where  V  is  the  bel  whose  direction  goes  from  v’  to  v  and  c{br )  is  the 
joint  cost  of  b 

(ii)  if  cctmp  <  cc{vr)  then 

a.  set  cc(v')  to  cctmp  and  dir(t;')  to  the  direction  from  v*  to  v; 

b.  if  v1  £  Q  then  insen  v'  in  Q  else  update  v '  in  Q ; 
endif 

endfor, 

endwkile; 

2.  starting  from  vc,  trace  recursively  the  next  vertex  pointed  to  by  the  current  vertex  using  the  direction 
information  in  dir  until  rs  is  reached,  and  return  the  bels  so  traced; 


end 

In  the  algorithms  above,  Q  is  a  bucket  represented  by  a  circular  vector  with  C  4*  1  positions  from  0  to  C  (see 
Figure  1).  Each  position  i,  i  =  0, . . . ,C,  has  associated  with  it  a  doubly  linked  list  which  contains  vertices  with  the 
same  cumulative  cost  value.  In  Step  3b(ii)c  of  algorithm  LWOF,  we  remove  all  vertices  v  from  Q  in  0(C)  time  since 
we  just  have  to  set  to  null  the  list  associated  with  each  position  t,  i  =  0, . . .  ,C,  in  Q.  An  index  io  is  used  to  indicate 
the  current  initial  position  in  Q  (see  Figure  1).  In  Step  la  of  algorithm  DSP,  a  vertex  v  in  Q  with  the  minimum 
cumulative  cost  cc(v)  is  removed  from  the  beginning  of  the  doubly  linked  list  at  position  to.  If  this  list  is  empty, 
to  is  incremented  until  the  next  position  in  which  a  non-empty  list  is  found.  Taking  the  worst  case,  this  operation 
has  a  computational  time  complexity  of  0(C) .  In  Step  lb(ii)b  of  algorithm  DSP,  a  vertex  v*  with  cumulative  cost 
cc(v')  is  inserted  in  Q  at  the  beginning  of  the  doubly  linked  list  at  position  [cc(z/)  mod  (C  +  1)].  This  operation  has 
a  computational  time  complexity  of  0(1).  The  Dijkstra’s  algorithm  guarantees  that  the  vertices  in  Q  will  be  always 
stored  in  the  increasing  order  of  cumulative  cost,  because  the  difference  between  the  maximum  and  the  minimum 


cumulative  costs  of  the  vertices  in  Q  is  always  less  than  or  equal  to  C.  In  the  same  step,  a  vertex  v'  in  Q  may  have 
its  cumulative  cost  updated,  meaning  that  we  have  found  a  new  path  from  va  to  v'  with  a  cost  less  than  the  current 
cost  cc(v ).  In  this  case,  we  have  to  remove  v‘  from  its  current  position  in  Q  and  insert  it  into  a  new  position  in  0 
This  process  is  done  with  a  computational  time  complexity  of  0(1). 


In  the  worn  case  algorithm  DSP  has  the  same  computational  time  complexity  0(m  +  nC)  as  in  the  Dial’s 
unplementation  of  Dykstra’s  algorithm,  where  m  is  the  number  of  bels  in  C,  n  is  the  number  of  vertices  in  C  and 

llZZ’TT  “n,  C(  )  Tf6*  *  beI  6-  0tber  s^rtest-path  algorithms  exist  with  computational  time 
leSVhp  °(m  +  7lC)  (ef’  °(m  +  nlogC),  0(m  +  nv^iZT),  and  0(m  log  log  C),  see20).  These  algorithms 
e  more  complex  data  structures  than  our  circular  queue  to  reduce  the  time  complexity  for  inserting  and  removing 

rr-*  OUr  mPlementation<  we  have  a  time  complexity  of  0(1)  for  inserting  and  updating  vertices  in  Q  In 
CaS?’  W6huVe  a  tune1compIexity  of  for  removing  a  vertex  from  Q  with  minimum  cumulative  cost,  as 
a  togamhnue  complexity  obtained  by  these  algorithms.  After  some  experimentation,  we  have  found  that 
e  number  of  increments  to  reach  the  next  non-empty  position  in  Q  is  usually  less  than  0.01  of  C.  Actually  even 
i  fe  Z  a,  ,g  Tabf  Typically>  we  have  used  4095  and  255  for  C  in  our  implementation  of  live  wire.  Therefore 
investigatldtethtr.  6  ^  Unpr0Vemeat “  Uve  wire  other  algorithms  is  really  significant.  This  should  be 


3.  EVALUATION 

as.sessed_the  g°°dness  of  a  segmentation  method  based  on  three  factors  -  precision,  accuracy,  and 
tte  r<2  I™  refers  to  the  repeatability  of  the  method  and  can  be  measured  by  evaluating  the  variations  in 
tn!th  EffiinrT' T tat*0n  because  of  subjective  operator  input.  Accuracy  refers  to  the  degree  of  agreement  with 
required  to ZZ  V°  *  y  °f  ^  meth°d  exPressed  35  some  Unction  of  the  total  user’s  time 

aXsfr  of  the 2  L  Serntf°n  Pl°CeSSu  d  °n  2’000  traCiDgS  “  a  P^^1^  ^PPhcation  and  statistical 
analysis  of  the  results,  we  have  shown  that  the  segmentations  of  the  2D  live-wire  method  in  general  agree  with 

those  of  manual  tracing  (accuracy)  and  that  the  live-wire  method  is  more  repeatable  (precision)  with  a  statistic^ 


significance  level  of  p  <  0.03,  and  1.5-2. 5  times  faster  (efficiency),  with  a  statistical  significance  level  of  p  <  0.02, 
than  the  manual  method.  In  this  section,  we  will  show  the  results  of  comparing  live  wire  and  live  wire  on  the  flv 
taking  into  account  the  efficiency  of  the  methods.  Since  the  delineation  of  the  contours  output  by  live  wire  on  the 
fly  is  exactely  in  the  same  way  as  in  live  wire,  the  accuracy  and  precision  of  the  former  will  be  identical  to  those  of 
the  latter,  and,  therefore,  they  need  not  be  assessed  again. 

In,6  we  have  introduced  a  feature  called  /8  in  live  wire  to  constrain  the  search  for  optimal  paths  in  the  current 
slice  to  an  annular  region  (shell)  of  width  W  centered  around  the  projection  onto  the  current  slice  of  the  contour(s) 
traced  in  the  previous  slice.  With  feature  /g,  live  wire  yields  very  fast  response  even  for  large  images.  Of  course,  we 
can  also  use  /8  to  further  improve  the  efficiency  of  live  wire  on  the  fly  on  large  images,  but  we  will  consider  in  this 
section  a  comparison  between  live  wire  with  /8  and  live  wire  on  the  fly.  Therefore,  our  experiments  will  take  into 
account  three  methods: 

•  LW:  live  wire  without  /s. 

•  LWF8:  live  wire  with  /g  using  W  =  60  pixels. 

•  LWOF:  live  wire  on  the  fly. 

For  our  experiments,  we  have  chosen  one  object  (the  talus  bone  of  the  human  foot)  in  one  of  our  ongoing 
applications,  the  kinematic  analysis  of  the  tarsal  joints  of  the  foot  based  on  MR  images.9”11  This  was  one  of  the 
objects  used  in  the  past  to  evaluate  the  previous  live  methods.6’16  We  created  a  set  of  67  2D  scenes  from  the  images 
within  our  database  as  follows.  The  images  (slices)  in  our  database  are  all  of  size  256  x  256  pixels.  We  chose  a  set, 
denoted  C256,  of  30  slices  from  this  set  pertaining  to  the  data  set  of  one  subject.  By  bilinear  interpolation  of  each 
of  these  slices,  we  created  another  set.  denoted  Ci2s,  of  30  128  x  128  slices.  Analogously,  we  created  a  set  C512  of 
five  512  x  512  slices  and  a  set  C1024  of  two  1024  x  1024  slices  from  the  original  256  x*256  slices.  The  reason  for 
choosing  a  fewer  number  of  scenes  of  size  512  x  512  and  1024  x  1024  is  that  the  response  time  of  LW  in  these  scenes  is 
prohibitively  slow.  One  operator  segmented  the  talus  in  each  of  these  scenes  using  each  of  the  two  methods  LW  and 
LWOF.  He  also  segmented  the  talus  in  C256  using  LWF8.  Our  evaluation  study  thus  consists  of  164  segmentation 
experiments  in  total.  More  experiments  involving  other  operators  are  currently  underway.  We  used  a  300  MHz 
Pentium  PC  for  these  experiments. 

We  denote  the  time  taken  to  complete  any  segmentation  experiment  e  by  Tt  (expressed  in  seconds).  Consider 
any  fixed  scene  type  t  €  {C128,  C256,  C512,  C1024}  and  method  m  €  {LW,  LWOF,  LWF8}.  We  define  the  time  taken 
Ttm  (in  seconds/slice)  for  segmenting  the  talus  in  a  2D  scene  of  type  t  using  method  m  to  be  the  average  of  all  times 

over  all  segmentation  experiments  e  involving  m  and  all  2D  scenes  of  type  t. 

We  have  done  three  types  of  timing  measurements.  The  first  type  measures  the  CPU  times  for  computing  the  live 
wire  segments  independent  of  other  supporting  processes  that  are  required  to  conduct  live  wire  segmentation.  This 
allows  us  to  assess  the  difference  in  speed  that  exists  purely  between  the  old  and  the  new  algorithms.  The  second 
type  measures  the  time  taken  by  the  user  to  segment  one  complete  contour  ignoring  the  time  for  other  processes  such 
as  displaying  the  slice  and  the  computation  of  the  cost  values  c{b)  for  all  bels.  The  third  type  includes  all  processes, 
and,  therefore,  gives  an  idea  of  the  comparative  user  time  required  for  overall  segmentation  for  the  different  methods 
in  an  actual  application.  We  note  here  that,  as  in  the  live-wire  method,6  training  is  required  only  once  for  an 
application  and  is  not  needed  on  a  per  study  basis.  This  is  typically  under  5  minutes  and  is  not  included  in  any  of 
the  time  measurements. 

Tables  1,  2  and  3  list  the  values  of  Ttm  for  all  possible  values  of  t  and  m,  for  the  three  types  of  measures, 
respectively.  Although  LWOF  can  find  optimum  paths  hundreds  of  times  faster  than  LW  (see  Table  1),  users 
cannot  react  with  the  same  speed  (see  Table  2).  Table  3  shows  that,  from  the  point  of  view  of  actual  segmentation, 
LW'  OF  is  about  1.5  to  33  times  faster  than  LWr  for  images  from  128  x  128  pixels  to  1024  x  1024  pixels.  Even 
constraining  optimum  path  computation  into  an  anular  region  of  width  equal  to  60  pixels  (i.e.,  method  LWrF$),  live 
wire  on  the  fly  is  about  2.3  times  faster.  Note  that,  the  advantage  of  live  wire  on  the  fly  over  live  wire  increases  with 
the  size  of  the  image  and  with  the  number  of  points  required  per  boundary.  In  our  experiments,  the  2D  boundaries 
of  the  talus  require  2  to  5  points  in  both  live  wire  strategies. 
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Table  1.  Segmentation  times 
of  measurement  that  indicates 


T  in  seconds /slice  for  all  possible  values  of  t  and  m.  This  table  hsts  the  first  type 
S^Haken  by  ,h,  shonest-path  algorithms  only  independent  of  other  processes. 
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T&ble  2.  Segmentation  times  Ttm 
type  of  measurement  that  indicates 
for  other  processes. 


in  seconds/slice  for  all  possible  values  of  t  and  m.  This  table  lists  the  second 
the  time  taken  by  the  user  to  segment  one  complete  contour  ignoring  the  tune 


Ci28 

C256 

C512  Cl024 

LW  ||  12 

24 

120  990 

LWOF  8 

8 

12  30  | 

LWF8  |  - 

|  18 

4.  CONCLUDING  REMARKS 

We  have  presented  a  new  user-steered  image  ^^^^^^f^re^^iy^ubSied  live  wbe  fraiSrork,6  but 
object  boundaries  in  a  slice-by-shce  fashion.  The  me  /soeed.  Based  on  164  segmentation  experiments 

utilizes  a  substantially  faster  shortest-path  gon  method  is  about  1  5  to  33  times  faster  than  live 

from  an  actual  medical  application,  we  have  shown  times  faster.  Other 

wire  for  actual  segmentation,  although  the  pure  computational  part  alone  is  over  a  nunor 

experiments  involving  multiple  operators  are  being  done.  ,  , 

A  drawback  of  live  ™  3  the  computation  tone  for  «  ^32 
on  the  boondary  to  all  other  points  in  the  image.  e  efficiency  considerably 

^activity ;of  the  method iln P^  — n^bl.Wsh^  ^  ^  ^  ^  aj  ^ 

» 5S? daring  the  segmentation  process.  Thus,  “he  «y  «"*“  “d 
Sl^e-eL  segments  in  real  tone,  even  for  very  large  images,  even  on  low-powered  comparers. 
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Abstract 

Studies  reported  in  the  literature  indicate  that  breast  cancer  risk  is  associated  with  mammographic 
densities.  An  objective,  repeatable,  and  a  quantitative  measure  of  risk  derived  from  mammographic 
densities  will  be  of  considerable  use  in  recommending  alternative  screening  paradigms  and/or  preventive 
measures.  However,  image  processing  efforts  toward  this  goal  seem  to  be  sparse  in  the  literature,  and 
automatic  and  efficient  methods  do  not  seem  to  exist.  In  this  paper,  we  describe  and  validate  an  auto¬ 
matic  and  reproducible  method  to  segment  dense  tissue  regions  from  fat  within  breasts  from  digitized 
mammograms  using  scale-based  fuzzy  connectivity  methods.  Different  measures  for  characterizing  mam¬ 
mographic  density  are  computed  from  the  segmented  regions  and  their  robustness  in  terms  of  their  linear 
correlation  across  two  different  projections  (CC  and  MLO)  are  studied.  The  accuracy  of  the  method  is 
studied  by  computing  the  area  of  mismatch  of  segmented  dense  regions  using  the  proposed  method  and 
using  manual  outlining.  A  comparison  between  the  mammographic  density  parameter  taking  into  account 
the  original  intensities  and  that  just  considering  the  segmented  area  is  demonstrated  which  indicates  that 
the  former  may  have  some  advantages  over  the  latter. 

Keywords:  Image  analysis,  mammograms,  glandular  tissue,  image  segmentation,  fuzzy  connectedness. 
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I.  Introduction 

In  the  mid  1970s,  studies  by  John  Wolfe  [1],  [2]  suggested  that  an  association  existed 
between  mammographic  parenchymal  patterns  and  the  risk  of  developing  breast  cancer. 
Since  then  there  have  been  many  studies  looking  at  the  relationship  between  mammo¬ 
graphic  fibroglandular  density  (often  referred  to  as  just  density)  and  the  risk  of  developing 
breast  cancer.  Although  a  few  studies  reported  no  association  of  density  with  increased 
risk,  the  majority  of  studies  have  found  an  association  between  parenchymal  patterns  and 
breast  cancer  risk.  A  recent  meta-analysis  [3]  of  all  studies  confirms  that  subjects  with 
mammographic  densities  have  an  increased  risk  of  breast  cancer  relative  to  those  without 
densities.  The  risk  increases  with  the  density  of  the  breast  [4].  Women  with  dense  breasts 
are  known  to  have  a  four  to  six  fold  increase  in  breast  cancer  risk  [1],  [5],  [6].  Cancers 
are  detected  at  later  stages  in  dense  breasts  and  mammographers  recognize  that  their 
diagnostic  accuracy  is  lower  in  such  women. 

The  Wolfe  classification  was  proposed  many  years  ago  to  identify  groups  of  women  at 
high  risk  for  breast  cancer  [5].  This  scheme  was  widely  used  for  many  years,  but  has 
fallen  into  disuse  because  of  several  limitations.  For  example,  inter-observer  variability  is 
a  problem  when  the  radiologists’  subjective  assessment  is  used  to  classify  the  amount  of 
density  present  [7] .  Secondly,  the  magnitude  of  the  increased  risk  has  varied  widely  in  the 
published  studies  [3].  Thirdly,  identification  of  this  risk  factor  for  a  given  woman  has  not 
altered  screening  recommendations  [6],  [8].  Computer-assisted  analysis  of  mammographic 
density  would  provide  an  objective,  quantitative  measure  of  cancer  risk  factor.  This  mea¬ 
sure  will  be  useful  in  total  risk  analysis  in  several  ways.  First,  such  risk  analysis  could 
influence  the  choice  of  alternative  screening  paradigms  such  as  intervals  between  mammo¬ 
grams  or  use  of  other  modalities  such  as  MRI.  Second,  this  measure  could  be  useful  in 
selecting  a  group  of  women  for  whom  the  risk-benefit  ratio  of  a  potentially  toxic  preven¬ 
tive  measure,  such  as  tamoxifin,  would  be  favorable  [16],  [17].  Third,  this  measure  could 
be  used  to  signal  the  need  for  more  careful  interpretation  of  a  subset  of  mammograms. 
For  example,  double-reading  might  be  indicated  for  mammograms  above  a  certain  level  of 
density. 

Image  processing  efforts  toward  this  goal  seem  to  be  sparse  in  the  literature,  and  au- 
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tomatic  and  efficient  methods  for  generating  this  measure  do  not  seem  to  exist.  Boyd 
et.  al.  [6]  studied  the  relation  between  mammographic  densities  and  breast  cancer  risk 
using  both  radiologist  classification  and  computer-assisted  density  measurement.  The 
computer-assisted  measurement  was  based  on  interactive  density  thresholding  using  two 
user-selected  thresholds.  They  observed  statistically  significant  increases  in  breast  cancer 
risk  associated  with  increasing  mammographic  density  in  both  methods.  Boone  et.  al. 
[18]  developed  and  evaluated  a  computerized  method  of  calculating  a  breast  density  index 
and  compared  this  index  with  breast  density  index  ranking  provided  by  mammographers. 
Byng  et.  al.  [19]  made  a  quantitative  symmetry  analysis  between  mammograms  of  differ¬ 
ent  breasts  of  the  same  patient  and  between  mammograms  at  different  projections  of  the 
same  breast  via  subjective  classification,  interactive  thresholding,  regional  skewness  mea¬ 
surement,  and  texture  analysis.  Ursin  et.  al.  [20]  studied  the  change  in  mammographic 
densities  in  women  participating  in  a  trial  of  a  gonadotropin-releasing  hormone  agonist 
(GnRHA)-based  regimen  for  breast  cancer  prevention  using  simultaneous  evaluation,  ex¬ 
pert  outlining,  and  non-expert  computer-based  thresholding  methods.  They  observed  that 
all  three  methods  yielded  statistically  significant  reduction  in  densities  from  baseline  to  the 
12-month  follow-up  mammogram  in  women  on  the  contraceptive  regimen.  They  found  a 
high  correlation  between  computer-based  results  and  the  results  from  the  expert  outlining 
method.  Huo  et.  al.  [21]  studied  the  ability  of  computer-extracted  features,  computed 
over  a  region  of  interest  selected  from  the  central  breast  region,  along  with  age,  to  identify 
women  at  risk.  They  found  that  a  computerized  characterization  of  parenchymal  patterns 
may  be  associated  with  breast  cancer  risk.  An  automatic  method  for  segmenting  the 
parenchymal  region  of  a  mammogram  using  first  and  second  order  gray-level  histograms 
is  presented  in  [22].  Another  approach  is  presented  in  [23]  for  determining  the  volume  of 
non  fatty  tissues  in  mammograms.  Recently  a  computer-assisted  user-interactive  method 
to  quantify  mammographic  density  has  been  published  [6],  which  concluded  that  quanti¬ 
tative  classification  of  densities  allows  for  the  determination  of  more  specific  gradients  of 
risk  than  do  Wolfe’s  classifications. 

In  this  paper,  we  describe  and  validate  an  automatic  and  reproducible  method  to  quan¬ 
tify  mammographic  densities  and  study  the  accuracy  of  related  parameters.  In  Section  II, 
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a  brief  description  of  the  principles  of  the  scale-based  fuzzy  connectedness  method  which 
forms  the  core  of  the  proposed  method  is  presented.  In  section  III,  we  describe  how  dif¬ 
ferent  parameters  are  automatically  selected  for  applying  fuzzy  connectivity  on  different 
regions.  In  Section  IV,  we  discuss  the  results  and  validate  the  method  —  (1)  by  studying 
linear  correlations  of  different  area  and  density  related  parameters  obtained  from  a  set  of 
mammograms  across  two  projections  and  (2)  by  studying  the  accuracy  of  the  method  by 
computing  the  area  of  mismatch  between  the  dense  regions  estimated  by  the  new  method 
and  by  manual  outlining.  A  comparison  between  the  mammographic  density  parameter 
taking  into  account  the  original  intensities  and  that  just  considering  the  segmented  area 
is  demonstrated.  Finally,  we  state  our  conclusions  in  Section  V. 

II.  Scale-Based  Fuzzy  Connectedness  Principles 

The  concepts  described  here  are  applicable  to  n-dimensional  (fuzzy)  digital  spaces;  see 
[24],  [25]  for  details.  However,  since  our  application  deals  with  two-dimensional  (2D) 
images,  we  confine  ourselves  only  to  the  2D  case. 

Most  real  objects  have  a  heterogeneous  material  composition.  Further,  imaging  devices 
have  inherent  limitations  including  spatial,  parametric,  and  temporal  resolutions.  In  the 
acquired  images  of  objects,  these  introduce  inaccuracies  and  artifacts  such  as  noise,  blur¬ 
ring,  and  background  variation.  The  artifacts  together  with  material  heterogeneity  cause 
the  object  regions  to  exhibit  a  gradation  of  intensity  values  in  the  image.  Even  if  the 
physical  object  is  perfectly  homogeneous  and  is  made  of  exactly  one  material,  its  image 
will  exhibit  a  graded  composition  within  the  object  regions  due  to  artifacts.  In  spite  of 
the  graded  composition,  knowledgeable  human  observers  usually  do  not  have  difficulties  in 
perceiving  object  regions  as  an  integrated  whole.  That  is,  image  elements  in  these  regions 
seem  to  hang  together  to  form  the  object  regions  in  spite  of  their  gradation  of  values. 
These  two  notions  —  graded  composition  and  hanging  togetherness  —  must  be  handled 
properly  by  any  segmentation  method  for  effective,  robust  performance.  In  our  methods, 
they  are  addressed  by  a  fuzzy  relation  among  image  elements  called  fuzzy  connectedness 
[24],  [25].  Such  a  general,  sound,  theoretic  and  algorithmic  framework  for  segmentation 
greatly  facilitates  the  quick  development  of  new  segmentation  applications,  as  we  have 
demonstrated  for  fuzzy  connectedness  in  brain  image  analysis  [10],  [11],  [12],  [13],  MRA 
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[14],  and  craniofacial  soft  tissue  display  [15]. 

The  scale-based  method  is  briefly  outlined  below  to  the  extent  needed  to  follow  our 
breast  segmentation  approach.  The  full  details  of  its  theory  are  given  in  [25].  We  will 
be  dealing  with  two  object  regions  in  our  segmentation  method  as  described  in  Section 
III.  The  first  corresponds  to  the  background  region  in  the  mammographic  image  and  the 
second  corresponds  to  the  dense  region.  In  the  description  in  the  rest  of  this  section, 
“object  region”  refers  to  each  of  these  regions. 

Throughout  we  denote  the  digitized  mammographic  image,  referred  to  as  a  (2D)  scene, 
by  C  =  ( C,f ),  where  C  denotes  the  pixel  array,  and  /(c)  denotes  the  pixel  value  for  any 
pixel  c  e  C.  For  any  pixel  c  £  C,  we  think  of  c  as  a  pair  (ci,c2)  representing  the  two 
coordinates  of  the  center  of  c.  The  range  of  /  is  assumed  to  be  [ L ,  H),  where  L  and  if  are 
integers. 

We  define  a  fuzzy  relation  k,  called  fuzzy  affinity,  on  the  pixel  array  C.  This  is  intended 
to  be  a  local  relation  among  pixels  that  are  nearby.  The  strength  of  this  relation  between 
any  two  pixels  c  and  d  in  C,  denoted  by  fxK(c,d),  lies  in  [0,1].  It  consists  of  three  com¬ 
ponents:  a  fuzzy  adjacency  component  a,  a  component  ip  based  on  object  homogeneity, 
and  a  component  <p  based  on  object  features.  The  idea  is  that  when  c  and  d  are  more 
adjacent,  have  more  homogeneity  of  intensities,  and  are  both  very  close  to  an  expected 
object  feature  value,  then  c  and  d  have  high  affinity.  In  other  words,  they  hang  together 
locally  very  strongly,  a  depends  on  how  far  c  and  d  are.  ip  depends  on  how  similar  the 
intensity  values  (or  other  features)  of  the  pixels  in  a  neighborhood  around  c  are  to  those 
around  d.  (p  depends  on  how  close  the  intensities  (or  other  features)  of  the  pixels  around 
c  and  those  around  d  are  to  some  expected  values  of  the  intensities  (or  other  features)  for 
the  object  under  consideration.  We  denote  the  strengths  of  all  these  three  components 
(all  of  which  lie  in  [0,1])  by  /xa(c,  d),  /x^(c,  d)  and  /^(c,  d),  respectively.  We  describe  below 
the  functional  forms  utilized  for  these  components. 

Although  the  theory  permits  more  general  forms,  in  this  paper,  we  use  the  following 
functional  form  for  /iQ.  For  any  pixel  c,  /ia(c.  c)  =  1.  Further,  for  any  two  pixels  c,d, 
Ha(c,d)  =  1  if  c  and  d  differ  in  exactly  one  coordinate  by  1;  otherwise  /xQ(c,  d )  =  0.  The 
specification  of  both  and  requires  the  notion  of  “object  scale”  at  every  pixel  in  C.  The 
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idea  behind  this  notion  is  that  if  we  can  roughly  estimate  the  size  of  the  object  structure 
locally  at  every  pixel,  then  this  information  can  be  utilized  to  determine  a  neighborhood 
size  around  c  and  d  for  specifying  ^  and  n#  in  a  way  that  is  tuned  to  the  object  and  is 
independent  of  pixel-level  variation  due  to  noise.  The  object  scale  r(c)  in  C  at  any  pixel 
c  in  C  denotes  the  size  (radius)  of  the  largest  disc  centered  at  c  that  lies  entirely  in  the 
object  region  in  which  c  lies.  Paradoxically,  it  appears  that  computing  scale  requires  image 
segmentation.  It  is  possible,  however,  to  develop  algorithms  that  give  a  rough  estimate 
of  object  scale  at  every  pixel  based  on  measuring  intensity  homogeneity  discontinuities 
and  that  do  not  require  explicit  image  segmentation.  We  have  demonstrated  in  [25]  that 
this  estimation  is  sufficient  to  give  a  good  approximation  of  the  scale  and  to  make  the 
fuzzy-connectedness-based  segmentation  very  robust  to  noise  and  pixel-level  variations. 
For  now,  we  assume  that  r(c)  is  known  at  any  c  €  C.  In  determining  the  scale-based  fuzzy 
affinity  between  any  pixels  c,  d  £  C,  two  digital  discs,  centered  at  c  and  d,  denoted  Bcd(c) 
and  Bcd{d),  both  of  radius  min[r(c),r(d)],  defined  by 

Bcd(c)  =  {eeC  |  ||c  -  e||  <  min[r(c),  r(d)]},  (1) 

Bcd(d)  =  {eeC\\\d-  e\\  <  min[r(c),  r(d)]},  (2) 


where  ||  •  ||  denotes  the  Euclidean  distance,  are  utilized. 

For  defining  m,  consider  any  two  pixels  c,  d  £  C  such  that  yua(c,  d)  >  0.  Consider  any 
pixels  e  £  BCd{c)  and  e'  £  Bcd{d )  such  that  they  represent  the  corresponding  pixels  within 
Bcd(c )  and  Bcd(d)-,  that  is,  c  —  e  =  d  —  e'.  We  will  define  two  weighted  sums  D+(c,  d)  and 
D~(c,d )  of  the  differences  of  intensities  between  the  two  discs  as  follows.  Let 


Scd(e > e') 


/(e)  -  Z(e'), 

0, 

V 

(  Z(e')  -  /(e), 

0, 


if  /(e)  -  Z(e')  >  0, 

otherwise, 

if  /(e)  -  Z(e')  <  0, 
otherwise. 


(3) 

(4) 


Then 


Y  [1  ))]  ^0,min[r(c),r(d)](||c  e||),  (5) 

e  G  Bcd(c) 
e'  G  Bcd(d) 
s.t.  c  —  e  —  d  —  g! 

"Y,  [1  —  Go ,171^+3^  (^cd(e)  6  ))]  Go,min[r(c),r(rf)]dlc  —  e||)>  (6) 

e  G  Bcd(c) 
e'  G  Bc<i(d) 
s.t.  c  —  e  —  d  —  e' 

where  “s.t.”  denotes  “such  that” ;  Gm,cT  denotes  an  unnormalized  Gaussian  with  mean  m 
and  standard  deviation  cr;  and  are,  respectively,  the  expected  mean  and  standard 
deviation  of  intensity  differences  between  all  pairs  of  adjacent  pixels  within  the  object 
region;  ||c  —  e||  represents  the  distance  between  c  and  e.  We  will  describe  in  the  next 
section  how  these  parameters  are  estimated  for  breast  images. 

The  connection  of  the  above  equations  to  the  homogeneity-based  affinity  ^  is  as  follows. 
There  are  two  types  of  intensity  variations  surrounding  c  and  d  —  intra-  and  inter-object 
variations.  The  intra-object  component  is  generally  random,  and  therefore,  is  likely  to  be 
near  0  overall.  The  inter-object  component,  however,  has  a  direction.  It  either  increases  or 
decreases  along  the  direction  given  by  c-d,  and  is  likely  to  be  larger  than  the  intra-object 
variation.  It  is  reasonable,  therefore,  to  assume  that  the  smaller  of  D+(c,d)  and  D~(c,d) 
represents  the  intra-object  component  and  the  other  represents  the  combined  effect  of  the 
two  components.  (Note  that  when  the  values  of  5^d  (respectively,  S~d)  are  small,  D+(c,d ) 
(respectively,  D~(c,d))  also  becomes  small.)  If  there  is  a  slow  background  component  of 
variation,  within  the  small  neighborhood  considered,  this  component  is  unlikely  to  cause  a 
variation  comparable  to  the  inter-object  component.  This  strategy  leads  us  to  the  following 
functional  form  for  ^ 

\D*(c,d)-D-[c,d)\ 

SeeBcd(c)  Go,min[r(c),r(d)](||c  —  e||) 

Note  that  \D+(c,  d)  —  D~(c,  d)\  represents  the  degree  of  local  inhomogeneity  of  the  regions 
containing  c  and  d.  Its  value  is  low  when  both  c  and  d  are  inside  an  (homogeneous)  object 


dip  (g  d)  —  1 


D+(c,d)  = 


ZT(c,d)  = 
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region.  Its  value  is  high  when  c  and  d  are  in  the  vicinity  of  (or  across)  a  boundary.  The 
denominator  in  (7)  is  a  normalization  factor. 

For  completing  the  specification  of  / 1 the  values  of  three  parameters  r(c),  and 
cty  need  to  be  determined.  The  method  of  estimating  r(c)  is  independent  of  the  type  of 
object  region  to  be  segmented  and  will  be  described  later  in  this  section.  The  method  of 
estimating  nhp  and  is  specific  to  the  object  region  and  will  be  explained  in  Sections 
I II- A  and  III-B. 

For  defining  the  object-feature-based  affinity,  we  first  compute  scale-based  filtered 
intensity  value  at  c  that  takes  into  account  the  disc  BT(c)  defined  by 

Br(c)  =  {e£C\\\c-e\\<r(c)}.  (8) 


The  filtered  intensity  value  at  any  c  e  C  is  given  by 

f  /  \  _  Se£gr(c)  /(e)  ^0,r(c)([lc  ~  ell)  /g\ 

EeeBr(c)Go>r(c)(||c-e||)  •  W 

Depending  on  whether  the  object  of  interest  is  darker  or  lighter,  we  define  the  following 
function  utilized  in  defining  /i^. 

Object  is  darker: 

( 

1, 


Object  is  lighter: 


W#(c)  = 


W*{c)  = 


i, 


(/.(c)), 

if  fa{c)  <  rritj,, 

otherwise. 

(10) 

(/.(c)), 

if  fa{c)  < 

(11) 

otherwise. 

In  both  cases,  m $  and  a#  represent  the  expected  mean  and  standard  deviation  of  the 
intensities  in  the  object  region,  respectively.  Finally,  the  following  functional  form  is  used 
for  nf 

1,  if  c  =  d, 

min[W^,(c),  W^(d)],  otherwise. 

For  completing  the  specification  of  n<p,  the  values  of  three  parameters  r(c),  and  cr# 
need  to  be  determined.  The  method  of  estimating  r(c)  is  presented  below.  The  method 
of  estimating  m </,  and  is  specific  to  the  type  of  object  region  to  be  segmented  and  will 


AV(C> d )  =  { 


(12) 
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be  explained  in  Sections  III- A  and  III-B.  Finally,  the  three  components  of  affinity  are 
combined  as  follows  to  define  the  scale-based  fuzzy  affinity, 

/j/c(c,  d)  —  Ha(c,  d)  n<t>(c,  d) .  (13) 

We  now  describe  an  algorithm  for  estimating  object  scale  at  every  pixel.  For  a  disc 
Bk(c)  of  any  radius  k  (see  (8))  centered  at  c,  we  define  a  fraction,  FOk{c),  that  indicates 
the  fraction  of  the  disc  boundary  occupied  by  a  region  in  which  the  scene  is  sufficiently 
homogeneous  with  c,  by 

t?/~ )  (  \ (c)—  B*_i(c)  £o,m*+3fl>  (|/(c)  -  /(d) |)  , 

\Bk(c)  —  Bk~i(c)\  •  1  j 

Here  \Bk(c )  -  Bk-i(c)\  denotes  the  number  of  pixels  in  Bk(c)  -  Bk-i{c).  We  define  B0(c) 
to  be  simply  the  set  {c}.  The  algorithm  for  object  scale  estimation  ( OSE )  is  summarized 
below.  The  algorithm  iteratively  increases  the  disc  radius  k  by  1,  starting  from  1,  and 
checks  for  the  fraction  of  the  object  FOfe(c)  containing  c  that  is  contained  on  the  disc 
perimeter.  The  first  time  when  this  fraction  falls  below  a  pre-selected  threshold  ts,  we 
consider  that  the  disc  enters  into  an  object  region  different  from  that  to  which  c  belongs. 
Following  the  arguments  in  [25],  we  have  used  ts  —  0.85. 

Algorithm  OSE 

Input:  C,  c  e  a  fixed  threshold  ts. 

Output:  r(c). 
begin 

set  k  —  1; 

while  FOk(c)  >  ts  do 
set  k  to  k  +  1; 
endwhile; 
set  r(c)  to  k ; 
output  r(c); 

end 

The  notion  of  pixel  affinity  captures  the  local  hanging-togetherness  property  of  pixels. 
The  notion  of  fuzzy  connectedness  expands  this  into  a  global  phenomenon  as  follows. 
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Consider  any  two  pixels  c  and  d  (not  necessarily  nearby)  in  C.  (Note  that  when  c  and 
d  are  far  apart,  their  affinity  is  0.)  Consider  any  path  (i.e.,  a  sequence  of  nearby  pixels) 
starting  from  c  and  ending  in  d.  We  define  a  “strength  of  connectedness”  of  this  path 
as  simply  the  smallest  affinity  (weakest  link)  along  the  path  between  a  pair  of  successive 
pixels  in  the  path.  Fuzzy  connectedness  is  a  global  fuzzy  relation,  denoted  K,  on  C.  The 
strength  of  this  relation  between  c  and  d  (not  necessarily  nearby),  denoted  pK(c,  d),  is  the 
largest  of  the  strength  of  connectedness  of  all  possible  paths  between  c  and  d.  A  scale- 
based  fuzzy  connected  object  of  C  of  strength  9,  for  any  9  €  [0, 1],  that  contains  a  specified 
pixel  o  in  C  is  a  subset  O  of  pixels  of  C.  O  is  such  that,  for  any  two  pixels  c,  d  in  O, 
px(c,  d)  >  9,  and  for  any  pixel  e  not  in  O ,  e)  <  9.  Given  C,  o,  9  and  a  scale-based 
affinity  relation  k,  finding  O  requires  the  computation  of  the  strength  of  connectedness 
of  literally  all  possible  paths  between  each  pair  of  pixels  in  the  set  of  all  possible  pairs  of 
pixels  in  C.  However,  the  theory  leads  to  practically  viable  algorithms  [24],  [25]  of  far  less 
complexity  that  are  based  on  dynamic  programming.  We  make  use  of  these  algorithms  in 
our  application. 

For  any  scene  C  =  (C,f),  any  fuzzy  affinity  k,  any  pixel  o  in  C,  we  define  the  fuzzy 
connectivity  scene  of  C  with  respect  to  o  to  be  the  scene  Cko  =  (C,  Jko),  where  for  any 
ce  C,  fKo(c)  =  I-Ik(o,  c ).  That  is,  the  value  assigned  to  any  pixel  c  in  CKo  is  the  strength 
of  connectedness  of  c  and  o.  We  generalize  this  definition  from  a  single  pixel  o  to  a  set  of 
pixels  X  by  setting  fKx(c)  =  ma xx€X{^k(x,c)}.  That  is,  in  the  fuzzy  connectivity  scene 
Ckx  =  (C,  !kx )  of  C  with  respect  to  X,  any  pixel  c  is  assigned  a  value  fKx  (c)  that  is  the 
maximum  of  the  strength  of  connectivity  of  c  with  the  elements  of  X.  Connectivity  scenes 
are  what  are  output  by  the  fuzzy  connectedness  algorithms  [24],  [25].  Upon  thresholding 
them,  we  get  the  segmented  fuzzy  objects. 

III.  Density  Quantification 

Our  method  of  mammographic  density  quantification  consists  of  the  following  steps: 
(1)  segmentation  of  the  breast  region  from  background;  (2)  segmentation  of  fat  and  dense 
regions  within  the  breast;  (3)  estimation  of  the  parameters  representing  quantified  density. 
These  are  described  in  separate  subsections  below. 
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A.  Segmentation  of  Breast  from  Background 

At  the  very  beginning,  using  3DVIEWNIX  [26]  supported  live-wire  [27]  tools,  regions 
corresponding  to  pectoral  muscles  are  interactively  excluded  when  those  are  projected  on 
to  the  scene.  This  tool  takes  help  from  the  operator  in  recognizing  where  the  pectoral 
muscles  are  in  the  image  but  does  the  delineation  of  their  boundary  automatically.  In  this 
fashion,  subjectivity  is  minimized.  In  the  entire  process  of  density  quantification,  this  is  the 
only  step  requiring  operator  intervention,  if  pectoral  muscles  appear  in  the  mammographic 
projection.  Scale-based  fuzzy  connectivity  is  used  for  segmenting  the  breast  region  from 
the  background.  Our  approach  will  be  to  segment  the  background  region  rather  than  the 
breast  region.  To  do  this,  we  need  to  (1)  determine  the  values  of  the  parameters  m$,  cr$, 
77ty,  and  <7 ^  for  the  background;  and  (2)  specify  a  set  of  pixels  in  the  background  region. 
These  are  accomplished  automatically  as  described  below. 

In  this  study,  we  have  utilized  120  mammograms  from  60  patients,  each  in  two  pro¬ 
jections  —  MLO  and  CC.  Studying  all  the  120  mammograms,  we  found  that  intensity 
histograms  of  mammograms  always  contain  a  prominent  peak  at  low  intensities,  and  this 
mode  corresponds  to  the  background.  A  typical  histogram  is  shown  in  Figure  1.  The  first 


Fig.  1.  A  typical  mammographic  intensity  histogram. 

prominent  peak  in  the  histogram  is  detected  and  the  intensity  corresponding  to  this 
peak  is  considered  as  the  mean  background  intensity.  Observing  that  this  part  of  the  his¬ 
togram  is  roughly  symmetric  about  m#,  the  standard  deviation  of  background  intensities 
00  is  computed  as  the  root-mean-squared  distance  of  the  intensities  from  m#  as  follows. 
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Let  h(i)  represent  the  number  of  pixels  in  the  mammographic  scene  with  intensity  i  (i.e., 
h(i)  is  the  height  of  the  histogram  at  i).  Then,  is  determined  as  follows 


— 


~  ^)2 /l(i) 


N 


(15) 


In  this  application,  since  the  object  of  interest  (i.e.,  the  mammographic  background)  is 
darker,  we  use  the  functional  form  of  (10)  for  computing  object-feature-based  affinity. 
Instead  of  an  operator  painting  pixels  in  the  background  region  for  training,  the  set  of 
pixels  in  C  satisfying  L  <  f(c)  <  m0  +  3(7#  is  utilized  for  estimating  the  parameters 
and  (Tip  for  homogeneity-based  affinity  /fy.  and  cty  are  taken  to  be  simply  the  mean 
and  the  standard  deviation  of  intensity  differences  between  all  pairs  of  adjacent  pixels  in 
this  set. 

During  digitization,  all  mammograms  were  oriented  so  that  the  top-right  and  the  bottom- 
right  corners  always  lie  in  the  background.  (This  was  standardized  for  both  left  and  right 
breasts.)  These  two  corner  pixels  are  used  to  form  the  reference  set  X  for  scale-based 
fuzzy  connectedness  processing.  Figure  2(b)  shows  the  connectivity  scene  obtained  for 
the  mammogram  at  CC  projection  shown  in  Figure  2(a).  As  shown  in  Figure  2(b),  there 
is  very  good  contrast  between  the  background  and  the  breast  region  in  this  connectivity 
scene.  We  discard  connectivity  strengths  greater  than  half  the  maximum  strength  and 
keep  the  lower  half  as  the  zone  for  the  breast  region.  This  zone,  however,  often  includes 
high  noise  pixels  and  markers  in  the  background  often  used  during  mammography.  To 
eliminate  these  pixels,  the  leftmost  1-pixel  in  the  middle  row  in  the  thresholded  connec¬ 
tivity  scene  is  chosen  as  the  reference  pixel  and  the  hard  connected  component  containing 
this  pixel  is  found  as  the  breast  region.  Figure  2(c)  shows  the  hard  segmented  breast 
region  for  the  original  mammogram  of  Figure  2(a).  This  method  has  worked  correctly  and 
automatically  in  all  studies  we  have  analyzed  so  far. 


B.  Segmentation  of  Fat  and  Dense  Regions 

Our  strategy  here  is  to  segment  the  dense  region  as  a  set  of  fuzzy  connected  objects.  The 
segmentation  operation  is  confined  to  the  breast  region.  The  fat  region  thus  gets  defined 
indirectly  as  the  complement  of  the  dense  region  in  the  breast.  For  this  segmentation,  as 
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(<*)  (6)  (c)  (d) 


Fig.  2.  (a)  An  original  mammographic  scene  of  a  patient’s  breast  at  CC  projection,  (b)  Scale-based 

fuzzy  connectivity  scene  for  the  background,  (c)  Segmented  breast  region,  (d)  Scale  based  fuzzy 
connectivity  scene  for  the  dense  region. 

in  breast  segmentation,  we  need  to  specify  the  values  of  parameters  m m^,  and  a $ 
for  the  dense  region  as  well  as  a  few  pixels  as  the  starting  information  in  the  dense  region. 

Our  approach  to  computing  the  parameters  will  be  as  for  the  segmentation  of  the  breast 
region  —  to  determine  automatically  a  set  of  pixels  that  are  definitely  in  the  dense  region 
and  then  to  estimate  the  parameters  from  the  intensity  distribution  within  this  set  of  pixels. 
In  this  application,  since  the  object  of  interest  (i.e.,  the  dense  regions)  is  lighter,  we  use  the 
functional  form  of  (11)  for  computing  object-feature-based  affinity  /^.  To  determine  a  set 
of  pixels  in  the  dense  region,  the  largest  intensity  value  MAX  is  determined  by  ignoring 
the  upper  0.1  percentile  of  intensity  in  the  histogram  of  the  breast  region.  Similarly, 
the  smallest  intensity  value  MIN  is  determined  by  ignoring  the  lower  0.1  percentile  of 
intensity.  We  then  select  within  the  breast  the  set  of  pixels  having  intensity  not  less  than 
MIN  +  0.75(MAX  —  MIN)  for  estimating  the  parameters.  Specifically,  we  take  the  mean 
and  standard  deviation  of  the  intensities  of  the  pixels  in  this  set  as  the  values  of  m#  and 
(7$.  Further,  we  take  the  mean  and  standard  deviation  of  intensity  differences  between  all 
pairs  of  adjacent  pixels  in  this  set  as  the  values  of  m ^  and  Finally,  the  set  of  pixels 
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in  the  breast  region  with  intensity  greater  than  MIN  +  0.85 (MAX  —  MIN )  is  used  as 
the  set  of  reference  pixels.  Figure  2(d)  shows  the  scale- based  fuzzy  connectivity  scene 
obtained  for  the  dense  region  for  the  original  mammogram  in  Figure  2(a).  In  the  next 
section,  we  describe  an  automatic  threshold  selection  method  that  is  applied  to  the  fuzzy 
connectivity  scene  for  the  dense  region  for  a  segmentation  of  the  breast  region  into  dense 
and  fatty  regions. 

C.  Automatic  threshold  selection 

The  connectivity  scene  is  an  image  in  which  pixels  hanging  together  strongly  with  the 
reference  pixels  supposedly  all  belong  to  the  same  object  region  [24].  It  is  usually  easy  to 
segment  the  object  region  by  thresholding  the  connectivity  scene  even  if  the  original  scene 
is  not  amenable  for  thresholding.  This  property  of  fuzzy  connectedness  is  demonstrated 
in  [24]  and  a  quantitative  validation  is  presented  in  [25].  In  the  past,  in  other  applica¬ 
tions  utilizing  fuzzy  connectedness  [10],  [13],  [14],  [15],  we  have  used  fixed  thresholds  on 
connectivity  scenes.  In  this  application,  we  found  that  fixing  the  connectedness  threshold 
is  not  always  satisfactory,  although  each  connectivity  scene  is  still  far  more  amenable  to 
thresholding  than  the  original.  With  this  as  the  motivation,  we  developed  an  automatic 
threshold  selection  strategy  for  connectivity  scenes,  building  on  the  idea  of  homogeneity- 
based  affinity  described  in  the  previous  selection  but  as  applied  to  the  connectivity  scene. 

To  select  the  best  threshold,  the  method  minimizes  an  energy  function  computed  by 
considering  spatial  arrangements  of  pixel  intensities  within  each  region  and  across  regions. 
We  emphasize  that  the  processing  is  now  confined  to  the  breast  region.  The  basic  idea  is 
as  follows.  Every  threshold  divides  the  scene  into  two  regions.  A  second  order  statistic, 
threshold  energy,  of  local  disagreements  in  the  scene  stemming  from  this  partitioning 
is  estimated  and  is  used  as  a  criterion  for  optimizing  the  threshold.  Threshold  energy 
characterizes  the  goodness  (rather,  badness)  of  a  particular  threshold  and  is  defined  as 
follows.  Let  B  denote  the  set  of  pixels  in  the  segmented  breast  region.  We  define  two 
fuzzy  relations  p  and  p,  respectively  called  likeliness  of  belonging  to  the  same  object  and 
likeliness  of  belonging  to  different  objects ,  on  the  pixels  in  B.  The  strengths  of  both  these 
relations  between  any  two  pixels  c  and  d  in  B  depend  on  (1)  how  far  c  and  d  are;  and 
on  (2)  how  similar  the  intensity  values  (or  other  features)  of  the  pixels  in  the  circular 
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neighborhood  around  c  are  to  those  around  d.  As  discussed  in  Section  II,  the  size  of  the 
circular  neighborhoods  around  c  and  d  depends  on  the  object  scales  at  c  and  d.  In  fact 
the  criterion  in  (2)  is  the  measure  of  homogeneity-based  affinity  between  c  and  d.  Two 
controlling  parameters  (indicating  expected  object  homogeneity)  are  required  to  calculate 
the  value  of  m  as  described  earlier.  These  two  parameters  are  estimated  as  the  mean  and 
the  standard  deviation  of  intensity  differences  of  all  pairs  of  adjacent  pixels  in  the  region 
with  pixel  intensities  (connectedness  values)  falling  in  the  upper  half  of  the  histogram  of 
the  connectivity  scene.  The  strength  of  the  fuzzy  relation  “likeliness  of  belonging  to  the 
same  object”  between  two  pixels  c,d  e  B,  denoted  fj,p(c,  d),  is  then  defined  as  follows. 


E 


AtP(c,  d)  =  /ia(c,  d) 


a,  b  £  B  s.t. 
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The  strength  of  the  fuzzy  relation  “likeliness  of  belonging  to  different  objects”  between 
two  pixels  c,d  e  B,  denoted  /ip(c, d),  is  defined  as  follows. 


d)  =  na(c,  d) 


1  <  MM) _ 

Eo,6eB  i  &) 
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Let  /p(c,  d,  t )  denote  a  predicate  that  takes  a  value  1  when  the  pixels  c,  d  belong  to  the  same 
object  at  the  threshold  t  and  0  otherwise.  Then  the  threshold  energy  E(t )  is  determined 
as  follows. 

E{t)  =  /P(c>rf^)Mp(c,d)  +  (1  ~  fp(c,d,t))fj,p(c,d).  (18) 

c,d€B 

In  words,  E(t)  expresses  the  level  of  concordance  between  the  two  regions  resulting  by 
applying  the  threshold  t  to  the  connectivity  scene.  Finally,  the  threshold  for  which  E(t) 
is  minimum  (indicating  minimum  concordance  or  maximum  discordance  between  the  two 
regions)  is  selected  as  the  optimum  threshold.  For  the  fuzzy  connectivity  scene  of  Figure 
2(d),  the  distribution  of  E(t)  is  shown  in  Figure  3(a),  while  Figure  3(b)  shows  the  location 
of  the  optimum  threshold  on  the  histogram  of  the  connectivity  scene  of  Figure  2(d).  The 
segmented  binary  scene  is  shown  in  Figure  3(c). 


Fig.  3.  (a),  (b)  Threshold  energy  and  connectivity  strength  distributions  for  the  connectivity  scene  of 

Figure  2(d).  (c)  Segmented  dense  region  using  automatic  thresholding. 

D.  Density  Quantification 

From  the  original  scene  and  the  segmented  fat  and  dense  regions,  the  following  param¬ 
eters  are  computed. 

TG:  Total  density  within  the  breast  region  defined  as  the  sum  of  intensities  of  pixels  in 
the  segmented  dense  region. 

TF:  Total  fat  within  the  breast  region  defined  as  the  sum  of  intensities  of  pixels  in  the 
segmented  fat  region. 

AB:  Total  area  of  breast  defined  as  the  number  of  pixels  in  the  segmented  breast  region. 
AG:  Total  area  of  density  within  the  breast  region  defined  as  the  number  of  pixels  in  the 
segmented  dense  region. 

AF:  Total  area  of  fat  within  the  breast  region  defined  as  the  number  of  pixels  in  the 
segmented  fat  region. 

AvF:  Average  pixel  intensity  within  the  fat  area  defined  by  TF/AF. 
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The  following  parameters,  some  of  which  are  derived  from  the  above,  which  may  be 
more  meaningful,  are  actually  used  in  our  testing:  TG,  TG/TF,  TG/AvF,  TG/AB,  AG, 
AG/AF,  and  AG/AB.  Linear  correlations  of  each  of  these  parameters  across  two  different 
projections  (CC  and  MLO)  were  tested  for  all  60  studies. 

IV.  Results  and  Discussion 

The  method  has  been  tested  in  two  ways  —  (1)  by  studying  correlation  of  each  density 
parameter  computed  from  patients’  mammograms  at  two  different  projections,  and  (2) 
by  studying  the  accuracy  of  the  method  in  terms  of  the  area  of  mismatch  between  dense 
regions  segmented  by  the  new  method  and  by  manual  outlining. 

A.  Patient  Mammograms 

The  method  has  been  tested  on  60  studies  selected  from  our  database.  Each  study  had 
two  mammographic  projections  —  CC  and  MLO.  These  mammograms  were  digitized  on  a 
Lumisys  scanner  at  a  resolution  of  100  microns.  The  population  includes  30  normal  studies 
as  well  as  14  studies  with  benign  and  16  with  malignant  masses  and  calcifications.  Except 
for  the  exclusion  of  pectoral  muscles  in  some  cases,  the  entire  method  worked  automati¬ 
cally  on  all  mammograms  wherein  all  parameters  required  by  the  algorithms  were  selected 
automatically.  An  additional  54  mammograms  were  processed  for  a  different  project  — 
to  assess  the  effect  of  hormone  therapy  on  breast  density.  The  algorithms  produced  vi¬ 
sually  acceptable  segmentations  in  all  174  cases.  Figure  4  demonstrates  the  results  of 
application  of  the  proposed  automatic  method  on  several  mammograms.  The  linear  corre¬ 
lation  coefficients  for  the  parameters  TG,  TG/TF,  TG/AvF,  TG/AB,  AG,  AG/AF,  and 
AG/AB  derived  from  the  two  sets  of  projection  images  were  0.967,  0.902,  0.951,  0.944, 
0.959,  0.915,  and  0.941,  respectively.  The  scatter  plots  and  the  R- values  of  different  pa¬ 
rameters  across  the  two  different  projections  over  60  pairs  of  mammograms  are  shown  in 
Figures  5(a)-(g).  In  all  these  figures,  the  horizontal  axis  represents  logarithmic  value  of 
the  estimated  parameter  for  each  mammogram  at  the  CC  projection  while  the  vertical 
axis  represents  the  same  for  the  matching  mammogram  at  MLO  projection.  Although  the 
scatter  plots  display  logarithmic  values,  the  R- values  of  linear  correlation  were  computed 
on  actual  values  of  the  parameters.  It  may  be  pointed  out  the  the  actual  scatter  in  the 


Fig.  4.  Results  of  application  of  the  proposed  density  segmentation  method  on  several  mammograms  at 
CC  and  MLO  projections.  In  each  set,  the  original  scene,  the  connectivity  scene  for  the  dense  region 
and  the  segmented  dense  region  are  shown. 

log-log  graphs  appear  less.  However,  the  log-log  graphs  were  used  to  present  compact 
displays  for  large  data  ranges.  The  ranges  of  gradients  and  y-intercepts  of  the  trend  lines 
were  0.998  to  1.0095  and  -0.2893  to  -0.0704,  respectively. 

The  high  value  of  correlation  coefficients  indicates  that  our  method  of  measurement  is 
highly  consistent  between  the  two  projection  images  of  the  same  patient.  The  highest 
correlation  is  obtained  for  TG.  Generally  the  parameters  that  use  area  measurements 
yielded  lower  correlations.  The  argument  behind  this  may  be  that  unless  the  3D  shape 
of  the  actual  dense  region  in  the  breast  is  approximately  spherical,  the  shapes  of  its 
CC  and  MLO  projections  may  be  quite  different  from  each  other.  This  may  yield  very 
different  area  measures  (AG,  AF)  although  the  total  density  may  still  be  the  same.  To 


(a) 


(h) 


(c) 


Fig.  5.  (a)-(g)  Scatter  plots  and  R- values  of  correlation  studies  for  different  quantitative  parameters 

computed  from  mammograms  at  CC  and  MLO  projections.  In  each  scatter  plot,  the  horizontal  axis 
indicates  the  logarithmic  value  of  the  parameter  at  CC  projection  while  the  vertical  axis  indicates  the 
same  at  MLO  projection. 

verify  this  hypothesis,  we  selected  among  the  60  pairs  of  studies  a  subset  of  20  pairs 
in  which  the  shapes  of  projections  of  the  same  breast  in  CC  and  MLO  appeared  quite 
different.  For  this  subset,  we  then  computed  the  correlation  coefficients.  The  coefficients 
for  TG  and  AG  for  this  subset  were  0.898  and  0.68,  respectively.  The  scatter  plot  for 
this  experiment  is  shown  in  Figures  6(a),(b).  To  estimate  statistical  significance  of  the 
difference  of  the  correlation  coefficients  for  these  two  populations,  we  used  the  Fisher’s 
Z-transformation  test  [30].  The  p- value  was  0.036245  showing  that  the  difference  in  the 
correlation  coefficients  is  statistically  significant.  However,  the  actual  data  were  skewed 
and  after  removing  the  stray  observations  the  p- value  rose  to  0.361511  indicating  that  the 
difference  was  not  significant.  The  primary  reason  behind  this  disagreement  may  be  the 
high  nonlinearity  between  the  total  length  in  3D  of  the  tissue  intercepted  by  a  beam  of 
X-ray  and  the  outcoming  energy.  If  this  nonlinearity  is  corrected  for  with  a  knowledge  of 
the  length  of  interception,  then  we  believe  that  TG  and  related  parameters  will  correspond 
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more  accurately  to  actual  total  density  than  the  area-based  parameters. 


(«)  (A) 


Fig.  6.  (a)-(b)  Scatter  plots  and  R-values  of  correlation  studies  of  TG  and  AG  for  a  set  of  20  studies 

with  quite  different  shape  at  different  projection.  In  each  scatter  plot,  the  horizontal  axis  indicates 
the  logarithmic  value  of  the  parameter  at  CC  projection  while  the  vertical  axis  indicates  the  same  at 
MLO  projection. 

B.  Comparison  with  Manual  Outlining 

A  major  difficulty  in  validating  a  density  quantification  method  is  how  to  generate  the 
truth.  Creation  of  a  realistic  physical  breast  phantom  with  known  volume  of  dense  tissue 
with  realistic  shape  and  distribution  is  a  research  topic  of  its  own.  In  this  paper,  we  have 
used  manual  outlining  of  dense  regions  by  an  expert  mammographer  to  provide  a  surrogate 
of  this  truth.  The  outlining  was  performed  using  a  computerized  freehand  drawing  tool 
on  digitized  mammograms  as  supported  by  3DVIEWNIX  [26].  An  example  of  manual 
outlining  is  presented  in  Figure  7.  The  task  of  manual  outlining  is  difficult,  quite  ill-defined, 
and  has  its  own  limitations  in  terms  of  definition  and  inter  and  intra-operator  variability. 
Additionally,  we  observed  that  the  correlation  of  the  density  parameters  computed  from 
manual  outlining  is  lower  than  that  using  the  proposed  method.  15  pairs  of  mammograms 
were  randomly  selected  from  our  data  set  of  60  pairs  of  mammograms.  Dense  regions  in 
each  mammogram  were  manually  outlined  by  the  same  expert  and  the  parameters  TG 
and  AG  were  computed  over  the  delineated  regions.  The  linear  correlation  coefficients  for 
TG  and  AG  were  0.638  and  0.534,  respectively.  The  scatter  plots  and  the  R-values  of 
these  two  parameters  across  the  two  projections  over  15  pairs  of  mammograms  are  shown 


Fig.  7.  An  example  of  manual  outlining.  The  border  of  the  dense  region  is  shown  bright. 


(«)  (b) 


Fig.  8.  (a)-(g)  Scatter  plots  and  R-values  of  correlation  studies  for  the  parameters  TG  and  AG  computed 
from  mammograms  at  CC  and  MLO  projections.  In  each  scatter  plot,  the  horizontal  axis  indicates 
the  logarithmic  value  of  the  parameter  at  CC  projection  while  the  vertical  axis  indicates  the  same  at 
MLO  projection. 

in  Figures  8(a)  and  (b). 

While  the  correlation  of  different  density  parameters  computed  by  the  proposed  method 
is  demonstrated  in  the  previous  section,  the  purpose  of  the  experiment  described  here  is 
to  show  that  the  disagreement  of  the  results  produced  by  the  automatic  method  with 
those  of  manual  outlining  is  within  the  limitation  range  of  the  second  method  itself.  15 
mammograms  were  randomly  selected  from  our  data  set  of  120  mammograms.  Dense 
regions  in  each  of  these  mammograms  were  manually  outlined  by  the  same  expert  at  two 
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different  time  instants  (with  a  gap  of  2  days).  Also,  the  dense  regions  were  segmented  in 
each  mammogram  using  the  proposed  method.  For  each  mammogram,  two  percent  areas 
of  mismatch,  AM,  were  computed,  one  between  the  results  of  manual  segmentation  at 
different  time  instants  and  the  other  between  the  results  of  manual  segmentation  at  the 
first  time  instant  and  those  using  the  proposed  method.  If  X  and  Y  are  the  two  sets  of 
pixels  being  compared,  AM(X,Y)  was  defined  by  x  100.  The  95%  confidence 

interval  of  AM  for  manual  segmentation  at  different  time  instants  was  [0,  22.81]  and  that 
for  manual  and  the  automatic  method  was  [0,  18.14].  This  shows  that  the  disagreement 
of  delineation  using  the  proposed  method  with  manual  outlining  is  within  the  range  of 
variability  of  the  latter  method  itself. 

V.  Conclusion 

A  near  automatic  method  for  quantification  of  breast  density  from  digitized  mammo¬ 
grams  has  been  developed  and  tested  on  87  pairs  of  patient  mammograms.  This  method 
was  executed  automatically  except  for  the  exclusion  of  projected  pectoral  muscles.  The 
method  consists  of  the  following  steps:  separation  of  the  breast  from  the  background,  cre¬ 
ation  of  a  fuzzy  connectivity  scene  for  the  dense  region,  segmenting  this  connectivity  scene 
using  an  automatic  threshold  selection  method,  and  then  computing  various  parameters 
that  characterize  total  breast  density.  A  set  of  density  and  area  related  parameters  has 
been  proposed  and  their  precision  in  terms  of  their  linear  correlation  across  two  different 
projections  has  been  studied.  The  scale-based  fuzzy  connectivity  method  has  been  found 
to  be  very  robust  and  effective  in  segmenting  the  mammographic  images.  Amount  of 
density  is  considered  to  be  one  of  the  strongest  risk  factors  for  breast  cancer.  Automatic, 
repeatable,  and  consistent  breast  density  quantification  from  digitized  mammograms  is 
practical  using  the  proposed  method.  The  correctness  of  the  proposed  density  quantifi¬ 
cation  method  has  been  validated  in  two  ways  —  (1)  showing  high  R-values  of  linear 
correlation  between  the  two  projections  (CC,  MLO)  of  the  various  parameters  computed 
over  segmented  dense  and  fatty  regions  and  (2)  demonstrating  the  disagreement  between 
delineations  using  the  proposed  method  and  using  manual  outlining  to  be  within  the  range 
of  variability  of  the  second  method.  The  method  removes  the  subjectivity  inherent  in  inter¬ 
active  threshold  selection  techniques  currently  used.  The  ability  of  the  computed  density 
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parameters  in  evaluating  risk  is  currently  being  investigated  at  our  institution. 

A  comparison  between  the  mammographic  density  parameter  taking  into  account  the 
original  intensities  and  that  considering  just  the  segmented  area  has  been  carried  out.  The 
motivations  behind  this  consideration  was  that  unless  the  3D  shape  of  the  actual  dense 
region  in  the  breast  is  approximately  spherical,  the  shapes  of  its  CC  and  MLO  projections 
may  be  quite  different  from  each  other  while  intensity  based  parameters  would  better 
capture  the  thickness  and  volume  of  the  dense  region  using  intensity  information.  For  60 
pairs  of  mammograms  at  CC  and  MLO  projections,  generally  the  parameters  that  use 
area  measurements  yielded  lower  correlations.  To  verify  this  argument,  we  selected  among 
the  60  pairs  of  studies  a  subset  of  20  pairs  in  which  the  shapes  of  projections  of  the  same 
breast  in  CC  and  MLO  appeared  quite  different  and  difference  in  the  correlation  coeffi¬ 
cients  was  statistically  significant.  However,  when  the  skewness  of  data  was  removed  by 
removing  stray  observations,  there  was  no  significant  difference  in  correlation  coefficients. 
The  primary  reason  behind  this  disagreement  may  be  the  high  nonlinearity  between  the 
total  length  in  3D  of  the  tissue  intercepted  by  a  beam  of  X-ray  and  the  outcoming  en¬ 
ergy.  Another  hurdle  in  using  intensity  related  parameters  originates  from  the  fact  that 
the  values  of  these  parameters  do  not  indicate  the  actual  volume  of  dense  regions.  This  is 
because,  actual  intensities  in  mammograms  are  dependent  on  different  imaging  parameters 
such  as  the  energy  and  the  frequency  of  X-rays  used,  plate  thickness,  film  characteristics 
and  X-ray  attenuation  coefficients  of  different  tissues  of  different  patients.  It  will  be  useful 
in  the  future  to  resolve  this  problem  of  nonlinearity  and  to  somehow  normalize  the  total 
density  parameters  such  that  they  relate  to  the  physical  volume  of  dense  regions. 

Acknowledgments 

The  authors  are  grateful  to  Drs.  Larry  Toto,  H.  L.  Kundel  for  the  digitized  mammo¬ 
graphic  data  set.  The  work  of  the  authors  is  supported  by  grants  DAMD  179717271  from 
the  Department  of  Army  and  NS  37172  from  NIH. 

References 

[1]  J.N.  Wolfe,  “Breast  pattern  as  an  index  of  risk  for  developing  breast  cancer”,  AJRy  126,  pp.  1130-1139, 
1976. 

[2]  J.N.  Wolfe,  “Breast  parenchymal  patterns  and  their  changes  with  age”,  Radiology ,  121,  pp.  545-552,  1976. 


r 


25 

[3]  E.  Warner  et.  al.,  “The  risk  of  breast-cancer  associated  with  mammographic  parenchymal  patterns:  a  meta¬ 
analysis  of  the  published  literature  to  examine  the  effect  of  method  of  classification” ,  Cancer  Detection  and 
Prevention ,  16,  pp.  67-72,  1992. 

[4]  A.M.  Oza  and  N.F.  Boyd,  “Mammographic  parenchymal  patterns:  a  marker  of  breast  cancer  risk”,  Epi¬ 
demiologic  Reviews ,  15,  pp.  196-208,  1993. 

[5]  J.N.  Wolfe,  “Risk  for  breast  cancer  development  determined  by  mammographic  parenchymal  pattern” ,  Can¬ 
cer ,  37,  pp.  2486-92,  1976. 

[6]  N.F.  Boyd,  J.W.  Byng,  R.A.  Jong,  E.K.  Fishell,  L.E.  Little,  A.B.  Miller,  G.A.  Lockwood,  D.L.  Tritchler, 
and  M.J.  Yaffe,  “Quantitative  classification  of  mammographic  densities  and  breast  cancer  risk:  results  from 
the  Canadian  national  breast  screening  study”,  Journal  of  the  National  Cancer  Institute  87,  pp.  670-675, 
1995. 

[7]  M.  Moscowitz,  P.  Gartside,  and  C.  McLaughlin,  “Mammographic  patterns  as  markers  for  high-risk  benign 
breast  disease  and  incident  cancers”,  Radiology ,  134,  pp.  293-5,  1980. 

[8]  J.N.  Wolfe,  A.F.  Saftlas,  and  M.  Salane,  “Mammographic  parenchymal  patterns  and  quantitative  evaluation 
of  mammographic  densities:  a  case  control  study”,  American  Journal  of  Roentgenology,  148,  pp.  1087-92, 
1987. 

[9]  J.  K,  Udupa,  L.  Wei,  Y.  Miki,  and  R.  I.  Grossman,  “A  system  for  comprehensive  analysis  of  multiple  sclerosis 
lesion  load  based  on  MR  imagery”,  SPIE  Proceeding ,  3031,  pp.  610-618,  1997. 

[10]  J.  K.  Udupa,  L.  Wei,  S.  Samarasekera,  Y.  Miki,  M.  A.  van  Buchem,  and  R.  I.  Grossman,  “Multiple  sclerosis 
lesion  quantification  using  fuzzy  connectedness  principles”,  IEEE  TRans.  Medical  Imaging ,  16,  pp.  598-609, 
1997. 

[11]  Y.  Miki,  R.  I.  Grossman,  J.  K.  Udupa,  M.  A.  van  Buchem,  L.  Wei,  M.  D.  Philips,  U.  Patel,  J.  C.  McGown, 
and  D.  L.  Kolson,  “Differences  between  relapsing  remitting  and  chronic  progressive  multiple  sclerosis  as 
determined  with  quantitative  MR  imaging” ,  Radiology ,  210,  pp.  769-774,  1999. 

[12]  Y.  Miki,  R.  I.  Grossman,  J.  K.  Udupa,  L.  Wei,  M.  Polansky,  L.  J.  Mannon,  and  D.  L.  Kolson,  “Relapsing- 
remitting  multiple  sclerosis:  longitudinal  analysis  of  MR  images  —  lack  of  correlation  between  changes  in 
T2  lesion  volume  and  clinical  findings”,  Radiology ,  213,  pp.  395-399,  1999. 

[13]  A.  Kumar,  W.  Bilker,  J.K.  Udupa,  and  G.  Gottlieb,  “Late  onset  minor  and  major  early  evidence  for  common 
neuroanatomical  substrates  detected  by  using  MRI”,  proceedings  of  the  National  Academy  of  Science ,  95, 
pp.  7654-7658,  1998. 

[14]  B.  L.  Rice,  Jr.  and  J.  K.  Udupa,  “Clutter-free  volume  rendering  for  magnetic  resonance  angiography  using 
fuzzy  connectedness”,  International  Journal  of  Imaging  Systems  and  Technology ,  11,  pp.  62-70,  2000.- 

[15]  J.K.  Udupa,  J.  Tian,  D.C.  Hemmy,  and  P.  Tessier,  “A  pentium-based  craniofacial  3D  imaging  and  analysis 
system”,  Journal  of  Craniofacial  Surgery ,  8,  pp.  333-339,  1997. 

[16]  S.M.  Lippman,  T.L.  Bassford,  and  F.L.  Meyskens,  Jr.,  “A  quantitatively  scored  cancer-risk  assessment  tool: 
its  development  and  use”,  Journal  of  cancer  Education ,  7,  pp.  15-36,  1992. 

[17]  R.T.  Chlebowski  et.  al.,  “Breast  cancer  chemoprevention  tamoxifen:  current  issues  and  future  prospective”, 
Cancer  72,  pp.  1032-7  Supplement,  1993. 

[18]  J.  M  Boone,  K.  K.  Lindfors,  C.  S.  Beatty,  and  J.  A.  Seibert,  “A  breast  density  index  for  digital  mammograms 
based  on  radiologists’  ranking”,  Journal  of  Digital  Imaging ,  11,  pp.  101-115,  1998. 

[19]  J.  W.  Byng,  N.  F.  Boyd,  L.  Little,  G.  Lockwood,  E.  Fishell,  R.  A.  Jong,  and  M.  J.  Yaffe,  “Symmetry  of 


v- 


26 

projection  in  the  quantitative  analysis  of  mammographic  images” ,  European  Journal  of  Cancer  Prevention , 
5,  pp.  319-327,  1996. 

[20]  G.  Ursin,  A.  Melvin,  A.  Astrahan,  M.  Salane,  Y.  R.  Parisky,  J.  G.  Pearce,  J.  R.  Daniels,  M.  C.  Pike,  and  D. 
V.  Spicer,  “The  detection  of  changes  in  mammographic  densities”,  Cancer  Epidemiology,  Biomarkers  and 
Prevention ,  7  pp.  43-47,  1998. 

[21]  Z.  Huo,  M.  L.  Giger,  O.  L  Olopade,  and  S.  A.  Cummings,  “Computerized  analysis  of  parenchymal  patterns 
for  the  assessment  of  breast  cancer  risk”,  Supplement  to  Radiology,  RSNA ,  209(P),  pp.  354,  1998. 

[22]  J.  Suckling,  D.  R.  Lewis,  and  S.  G.  Blacker,  “Parenchymal  delineation  by  human  and  computer  observers”, 
in  Proceeding  of  the  2nd  International  Workshop  on  Digital  Mammography,  Eds.  A.G.  Gale  et.  al.,  York, 
England,  1994. 

[23]  R.  Highnam,  M.  Brady,  and  B.  Shepstone,  “A  representation  for  mammographic  image  processing”,  Medical 
Image  Analysis ,  1,  pp.  1-18,  1996. 

[24]  J.  K.  Udupa  and  S.  Samarasekera,  “Fuzzy  connectedness  and  object  definition:  theory,  algorithms,  and 
applications  in  image  segmentation”,  Graphical  Models  and  Image  Processing ,  58,  pp.  246-261,  1996. 

[25]  P.  K.  Saha,  J.  K.  Udupa,  and  D.  Odhner,  “Scale-based  fuzzy  connected  image  segmentation:  theory,  algo¬ 
rithms,  and  validation”,  Computer  Vision  Image  Understanding ,  77,  pp.  145-174,  2000. 

[26]  J.K  Udupa,  D.  Odhner,  S.  Samarasekera,  R.J.  Goncalves,  K.  Iyer,  K.  Venugopal,  and  S.  Furuie, 
“3DVIEWNIX:  a  open,  transportable,  multidimensional,  multimodality,  multiparametric  imaging  system, 
proceedings  of  SPIE ,  2164,  pp.  58-73,  1994. 

[27]  A.X.  Falcao,  J.K.  Udupa,  S.  Samarasekera,  and  S.  Sharma,  “User-steered  image  segmentation  paradigms: 
live  wire  and  live  lane”,  Graphical  Models  Image  Processing  60,  pp.  233-260,  1998. 

[28]  M.  Tabb  and  N.  Ahuja,  “Multiscale  image  segmentation  by  integrated  edge  and  region  detection”,  IEEE 
Trans .  Image  Processing ,  6,  pp.  642-655,  1997. 

[29]  T.  Lindeberg,  Scale- Space  Theory  in  Computer  Vision ,  Boston,  MA:  Kluwer,  1994. 

[30]  R.  F.  Woolson,  Statistical  methods  for  the  analysis  of  biomedical  data ,  New  York:  Wiley,  1987. 

[31]  T.S.  Curry,  J.E.  Dowdey,  and  R.  C.  Murry  Jr.,  Christensen’s  Introduction  to  the  Physics  of  Diagnostic 
Radiology ,  Lea  &;  Febiger,  Philadelphia,  1984. 


